function [y,s,means] = remmed(X,d,meanmed);
%REMMED Remove mean- or median-based seasonal component.
%   Y = REMMED(X,D) returns time series X with removed mean-based seasonal 
%   component of period D. E.g. REMMED(X,7) returns daily data without 
%   the mean week as described in [1], Section 2.4.2.
%   Y = REMMED(X,D,1) returns time series X with removed median-based 
%   seasonal component.  
%   [Y,S,MEANS] = REMMED(X,D) additionally returns the seasonal component 
%   S and the mean level MEANS of the seasonal component S.
%
%   Reference(s):
%   [1] R.Weron (2007) 'Modeling and Forecasting Electricity Loads and 
%   Prices: A Statistical Approach', Wiley, Chichester.   

%   Written by Rafal Weron (2005.08.07)
%   Copyright (c) 2005-2006 by Rafal Weron

if nargin<3,
   meanmed = 0;
end;

% Make a column vector
X = X(:);

% Make length X a multiple of d
N = length(X);
D = floor(N/d);
x = X(1:D*d); 

% Reshape data
rx = (reshape(x,d,D))';

if meanmed == 1,  
    mx = median(rx);    % meadian
else            
    mx = mean(rx);      % mean
end

% Seasonal component
means = mean(mx);
s = mx' - means;
S = repmat(s,D,1);

% Remove seasonal component
y = X - [S; S(1:N-D*d)];