clear all;
close all;

%%%
% Load data
%%%
imp = importdata('auto1.asc');
y  = imp.data(:,end);
x2 = imp.data(:,end-2);
x3 = imp.data(:,end-1);

%%%
% Calendar utilities
%%%
N = length(y);
cl=cal(1959,1,4,N);
stop=ical(1973,3,cl);
yr=tsdate(1973,1,4,N);
years =59:0.25:(59+(N-1)*.25);

%%%
% OLS estimation
%%%
X = [ones(N,1) x2 x3];
[N k]=size(X);
stop=stop-8;
reg = ols(y(1:stop),X(1:stop,:));
prt(reg)

%%%
% Plot data
%%%
yr=y(1:stop);
x2r=x2(1:stop);
x3r=x3(1:stop);
yr = (yr-mean(yr))/std(yr);
x2r = (x2r-mean(x2r))/std(x2r);
x3r = (x3r-mean(x3r))/std(x3r);
lgd=['Gas  ',
     'Inc  ',
     'Price'];
figure(1)
tsplot([yr x2r x3r],cl,1,stop,lgd);
figure(2)
pairs([years' y x2 x3],['Years';lgd])
figure(3)
gplotmatrix([years' y x2 x3],[],ones(N,1),[],[],[],[],'hist',{'Years';'Gas';'In';'Price'},{'Years';'Gas';'In';'Price'})

%%%
% Plot recursive residuals
%%%
N=stop+8;
cs = cusums(y(1:N),X(1:N,:));
sg = zeros(N,1);
for j=k+1:N
    reg1 = ols(y(1:j),X(1:j,:));
    e = reg1.resid;
    sg(j) = sqrt(e'*e/(j-k)); 
end
figure(4);
plot(years((k+1):N),cs.rres((k+1):N),'-r',years((k+1):N),2*sg((k+1):N),':b',years((k+1):N),-2*sg((k+1):N),':b')
axis([years(k+1) years(N) -0.06 0.08])

%%%
% Plot scaled F-test values for one-step forecast error
%%%
RSS = zeros(N,1);
Fstat = RSS;
betaj=zeros(N,k);
betaj_lo=betaj;
betaj_up=betaj;
reg1 = ols(y(1:(k+1)),X(1:(k+1),:));
e = reg1.resid;
RSS(k+1) = sqrt(e'*e);
for j=k+2:N
    regj = ols(y(1:j),X(1:j,:));
    e = regj.resid;
    betaj(j,:)=(regj.beta)';
    betaj_lo(j,:)=(regj.bint(:,1))';
    betaj_up(j,:)=(regj.bint(:,2))';
    RSS(j) = sqrt(e'*e); 
    Fstat(j)=2.1*(j-k-1)*(RSS(j)-RSS(j-1))/RSS(j-1)/finv(0.95,1,j-k-1);
end
figure(5);
plot(years((k+1):N),Fstat((k+1):N))
line([years(k+1) years(N)], [1 1],'LineStyle','--')
axis([years(k+1) years(N) 0 3.3])

%%%
% Recursively estimated coefficient
%%%
figure(6)
dlt=7;
h=1;
subplot(2,2,h)
plot(years((k+dlt):N),betaj((k+dlt):N,h),'-r',years((k+dlt):N),betaj_lo((k+dlt):N,h),...
    ':b',years((k+dlt):N),betaj_up((k+dlt):N,h),':b')
ylabel('Constant')
xlabel('Year')
axis tight
h=2;
subplot(2,2,h)
plot(years((k+dlt):N),betaj((k+dlt):N,h),'-r',years((k+dlt):N),betaj_lo((k+dlt):N,h),...
    ':b',years((k+dlt):N),betaj_up((k+dlt):N,h),':b')
ylabel('Price elasticity')
xlabel('Year')
axis tight
h=3;
subplot(2,2,h)
plot(years((k+dlt):N),betaj((k+dlt):N,h),'-r',years((k+dlt):N),betaj_lo((k+dlt):N,h),...
    ':b',years((k+dlt):N),betaj_up((k+dlt):N,h),':b')
ylabel('Income elasticity')
xlabel('Year')
axis tight

%%%
% CUSUM and CUSUMSQ tests of Brown,Durbin and Evans
%%%
figure(7)
wt = cs.rres;
wt(1:k)=0;
Wt=cumsum(wt);
Wt=Wt/sqrt(reg.sige);
plot(years(1:N),Wt,'-r')
line([years(k) years(N)],[ 0.948*sqrt(N-k)  3*0.948*sqrt(N-k)],'Color','b','LineStyle',':')
line([years(k) years(N)],[-0.948*sqrt(N-k) -3*0.948*sqrt(N-k)],'Color','b','LineStyle',':')
line([years(k) years(N)],[0 0],'Color','k','LineStyle',':')
title('CUSUM with 95% intervals');
legend('cusum','upper95','lower95','Location','NorthWest');
xlabel('Observations');

figure(8)
c0tab=[.475 .50855 .46702 .44641 .42174 .40045 .38294 .36697 .35277 .34022 .32894 .31869 ...
    .30935 .30081 .29296 .28570 .27897 .27270 .26685 .26137 .25622 .25136 .24679 .24245...
    .23835 .23445 .23074 .22721 .2238 .2206 .2175 .2146 .2117 .2090 .2064 .2039 .2014...
    .1991 .1968 .1946];
[coef,r,j]=nlinfit(10:40,c0tab(10:40),@c0fun,[2 2 .5]);
plot(years(1:N),cs.cusums,'-r')
c0=c0fun(coef,N);
line([years(k) years(N)],[c0 1+c0],'Color','b','LineStyle',':')
line([years(k) years(N)],[-c0 1-c0],'Color','b','LineStyle',':')
line([years(k) years(N)],[0 1],'Color','k','LineStyle',':')
title('CUSUMSQ with 95% intervals');
legend('cusumsq','upper95','lower95','mean value','Location','NorthWest');
xlabel('Observations');