addpath ../Datasets
clear; 
%close all; utl_mfig;
clf

% Load Keeling Curve
load CO2.mat
t=[1:length(co2)]';
y=co2(~isnan(co2));
mytime=mydate(~isnan(co2));
t=t(~isnan(co2));	
%std_n=2;
std_n=0.1;  % measurment error


N=length(t);

% Observation Error Covariance
Rnn = diag( ones(N,1)*(std_n).^2 );

% Set the matrix  y = a + b*t + c*t^2 + d*cos(2*pi*t/12) + e*sin(2*pi*t/12)
E = ones(N,5);
E(:,1)=1; E(:,2)=t; E(:,3)=t.^2; E(:,4) = cos(2*pi*t/12); E(:,5) = sin(2*pi*t/12); 

param=[1 2 3 4 5];
%param=[1 2 3];

E=E(:,param);


% Estimate xhat
CME=inv(E'*E);
xhat = (CME)*E'*y
% Uncertainty in model parameters
P = CME*E'*Rnn*E*CME;
xhat_uncer = diag(P)

% Define matrix H
H=E*CME*E';
yhat = E*xhat; % Equivalent to H*y
nhat = y - yhat;

J = length(t) - size(E,2)
Jhat=nhat'*nhat/(std_n^2)

disp([' Residuals are orthonogal to noise free estimate =', num2str(nhat'*yhat)]);
subplot(2,2,1)
plot(mytime,y,'k'); hold on
plot(mytime,yhat,'r');
legend 'OBSERVATIONS' 'LSQ ESTIMATE' 
datetick
title 'LSQ fit'
ylabel 'PPM'
xlabel 'YEAR'
subplot(2,2,3)
plot(mytime,nhat,'k'); hold on
%plot(mytime,nhat,'k.');
datetick
title 'RESIDUALS'
ylabel 'PPM'
xlabel 'YEAR'

subplot(2,2,2)
modparam=[1 2 3];
y_trend = E(:,modparam)*xhat(modparam);
plot(mytime,y_trend,'k'); hold on
title 'TREND'
datetick
ylabel 'PPM'
xlabel 'YEAR'

if length(param)>3
subplot(2,2,4)
modparam=[4 5];
y_seas = E(:,modparam)*xhat(modparam);
plot(mytime,y_seas,'k'); hold on
title 'SEASONAL'
datetick
ylabel 'PPM'
xlabel 'YEAR'
end