function [y,s,P] = remst(X,d,degree,glue);
%REMST Remove seasonal component and trend.
%   Y = REMST(X,D,DEGREE) returns time series Y with removed polynomial 
%   trend of degree DEGREE and seasonal components of period D from vector 
%   X. REMST uses the moving average technique (see [1], Section 2.4.3).
%   Y = REMST(X,D,DEGREE,1) glues the first and last D/2 values to the end 
%   and the beginning of X, respectively (recommended for short time 
%   series, i.e. of length 2-3 D).
%   [Y,S,P] = REMST(X,D,DEGREE) returns seasonal component S and 
%   polynomial coefficients P.
%
%   If DEGREE==-1 then only seasonal components are removed. 
%   If DEGREE==-2 then only seasonal components are removed and the minimum 
%   values of the original and deseasonalized series are aligned.
%   If DEGREE is a vector it is treated as a one-period (i.e. of length D) 
%   seasonal component and is subtracted from X.
%
%   Reference(s):
%   [1] R.Weron (2007) 'Modeling and Forecasting Electricity Loads and 
%   Prices: A Statistical Approach', Wiley, Chichester.   

%   Written by Joanna Nowicka-Zagrajek and Rafal Weron (2001.01.27, rev. 2006.09.23)
%   Copyright (c) 2001-2006 by Rafal Weron

if nargin<4,
    glue = 0;
end;

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

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

% Perform computations if no seasonal component was provided
if length(degree)<2,
    if glue == 1, % glue first and last d/2 values to the series
        if (d/2) == floor(d/2)
            q = d/2;
        else
            q = (d-1)/2;
        end;
        x = [X(1:D*d)' X(1:2*q+1)']';
    else % use the original series
        x = X(1:D*d);
    end;

    n = length(x);
    m = zeros(n,1);

    % Calculate mean
    if (d/2) == floor(d/2), % even
        q = d/2;
        for t = q+1:n-q
            m(t) = sum([0.5*x(t-q) x(t-q+1:t+q-1)' 0.5*x(t+q)])/d;
        end;
    else % odd
        q = (d-1)/2;
        for t = q+1:n-q
            m(t) = mean(x(t-q:t+q));
        end;
    end;
    x = x(q+1:n-q);
    m = m(q+1:n-q);

    % Subtract the mean
    xm = x - m;
    w = zeros(d,1);
    for k = 1:d
        w(k) = mean(xm(k:d:end));
    end;

    % Compute seasonal component for one period
    s = zeros(d,1);
    for k = 1:d
        s(mod(k+q-1,d)+1) = w(k) - mean(w);
    end;

else  % DEGREE is the one-period seasonal component  
    s = degree;
    s = s(:);
end % if (length(degree)<2)

% Compute seasonal component for the whole series
ss = zeros(N,1);
for i = 0:D-1,
    ss(i*d+1:i*d+d) = s;
end;
ss(d*D+1:N) = s(1:(N-D*d));

% Perform computations if no seasonal component was provided
if (length(degree)<2),
    if degree < 0, 
        y = X - ss;
        P = [];
    else
        T = 1:N;
        P = polyfit(T',X-ss,degree);

        y = zeros(N,1);
        for t = 1:N,
            y(t) = X(t) - ss(t) - dot(P(degree+1:-1:1),[1 cumprod(ones(1,degree)*t)]);
        end;
    end;
    if degree == -2,
        y = y + min(X) - min(y);
    end
else % DEGREE is the one-period seasonal component 
    y = X - ss;
    P = [];
end