%%  This modification is for the case when we first fit a trend 
%%  --here, quadratic,
%%  and then fit the seasonal. We also attempt to find autocorrelation

X = load('./data/mauna-loa.csv');

co2 = X(:,2);
co2 = co2(11:length(co2));   % throw out first 10 months 

n = length(co2);  % number of obs 

years = 50;

M = 6;

w = ones(2*M+1,1)/12;
w(1) = 1/24;
w(13) = 1/24;

m = co2*0;
for t = (M + 1):(n-M),
   m(t) = sum( co2((t-M):(t+M)).*w); 
  
end

plot(X)
  hold on
plot(m,'r'); 











% =========================================================================
% Part 1. Fitting a trend.
                    % fit a second-order trend 

X = (1:n)' - n/2;        % years where 1 = 1959, then centered to avoid computational difficulties 
G = [ones(n,1) X X.^2];

Y = co2;
betahat = (G' * G)^(-1) * G' * Y;      % the least squares solution
Ypred = G * betahat;


plot(Y) 
  hold on
  plot(Ypred,'r')
  hold off

resid = Y - Ypred;
s = sqrt( sum(resid.^2)/(n-3));   % this is not yet good because the seasonal component remains
 plot(resid);
 
% =========================================================================
% Part 2. Estimating seasonal fluctuations
  
 resmo = reshape(resid, 12, years) ;      % make it a rectangular array 
 seas = mean(resmo');                     % take average over each month
 plot(seas,'*')                   % seasonal profile
 
 seas50 = seas' *ones(1,50);     
 seas_rep = reshape(seas50,600,1);    % repeated 50 times
 plot(seas_rep) 

 Yclean = Y - seas_rep;     % data without seasonal fluctuation

 betahat2 = (G' * G)^(-1) * G' * Yclean;      % the least squares solution based on cleaned data
 Ypred2 = G * betahat2;

  res2 = Yclean - Ypred2;
  s2 = sqrt( sum(res2.^2)/(n-3));   % this one's better since we removed the seasonal component
  plot(res2);
  
         % this is not yet very satisfying because some parts of the trend remain.  
         
  
% =========================================================================
% Part 3. Modeling autocorrelation

 c = xcov(res2);      % xcov is "cross-covariance" between vectors, or autocovariance in terms of one vector
 plot(c)               % function xcorr is called "cross-correlation" by Matlab, but the name is misleading
 
 c = c/(s2^2*n);
 plot(c)             % this is what *we* call autocorrelation   
   title 'Autocorrelation function, lag 0 = 600';

  maxlag = 300; 
  lag = 0:maxlag;
  plot(lag,c(599:(599+maxlag))) 
   
                 % you can still see some artifacts:  fairly large
                 % autocorrelation at the lag 300 
 
 
% ==================== 
% ==============
% Part 4. Difference modeling


dY = res2(2:600) - res2(1:599);
plot(dY)
plot(xcov(dY))

c = xcov(dY);
plot(c(599:605))