%Macro 202a Homework #2 Solution Anton A Cheremukhin 19 October 2005
close all;  clear all;  clc;
%data
[time, gdp, cons, inv1, hours,inv2] = textread('datarbc.txt','%d %f %f %f %f %f');
T=size(cons,1);
gdphp=hpfilter(gdp,100);    y0=log(gdp./gdphp);    
conshp=hpfilter(cons,100);  c0=log(cons./conshp);
hourshp=hpfilter(hours,100);    l0=log(hours./hourshp);

%c0=log(cons);   X=[ones(T,1) time]; beta=inv(X'*X)*(X'*c0); c0=c0-X*beta;
%y0=log(gdp);   X=[ones(T,1) time]; beta=inv(X'*X)*(X'*y0); y0=y0-X*beta;
%l0=log(hours);   X=[ones(T,1) time]; beta=inv(X'*X)*(X'*l0); l0=l0-X*beta;

[capital] = textread('capital.txt','%f');
capitalhp=hpfilter(capital,100);    k0=log(capital./capitalhp);    
%k0=log(capital);   X=[ones(T,1) time]; beta=inv(X'*X)*(X'*k0); k0=k0-X*beta;

%predefined parameters
beta=0.96;      g=0.0018;        ro=3;
alpha=1/3;      fi=0.17;         lambda=0.65;    delta=0.1;
sigmay=0.0004;     sigmak=0.003;     sigmac=0.0001;    sigmal=0.001;
vol=0.0002;
%set of parameters
q0=[log(-1+1/beta), log(-1+1/g), log(ro), log(-1+1/alpha), log(-1+1/fi) ,log(-1+1/lambda), log(-1+1/delta), log(sigmay), log(sigmak), log(sigmac), log(sigmal), log(vol)];
y2=[y0,k0,c0,l0];

H0=0.0001*eye(12);

crit=0.00001;
nit=1000;
% find likelihood maximizer with csminwel
[fh,xh,gh,H,itct,fcount,retcodeh] = csminwel(@neolike, q0, H0, [], crit, nit, y2);

beta=1/(1+exp(xh(1)));
g=1/(1+exp(xh(2)));
ro=exp(xh(3));
alpha=1/(1+exp(xh(4)));
fi=1/(1+exp(xh(5)));
lambda=1/(1+exp(xh(6)));
delta=1/(1+exp(xh(7)));
sigmay=exp(xh(8));
sigmak=exp(xh(9))+0.00000000001;
sigmac=exp(xh(10));
sigmal=exp(xh(11));
vol=exp(xh(12));
%%%%%%%%%%%%%%%%%%%
%new parameters - definitions
teta=((((1+g)^(1-fi*(1-ro)))/beta)-(1-delta))/alpha;%%=cap/gdp
nu=1-(delta+g)/teta; %=cons/gdp
labor=1-1/(1+fi*(1-alpha)/(nu*(1-fi)));
%steady state values
cap=labor/(teta^(1/(1-alpha)));
gdp=teta*cap;
cons=nu*gdp;
%more parameter definitions
phi=alpha*beta*gdp/((1+g)*cap);
omega=fi*(1-ro)-1;
a1=labor/(1-labor);
pi=(1-fi)*(1-ro)*a1;
%matrices
D=[ 0,-phi,-omega,pi,0,phi,omega,-pi;
    -gdp,(1+g)*cap,cons,0,0,0,0,0;
    1,0,0,-(1-alpha),-1,0,0,0;
    1,0,-1,-(1+a1),0,0,0,0;
    0,0,0,0,1,0,0,0;
    1,0,0,0,0,0,0,0;
    0,1,0,0,0,0,0,0;
    0,0,1,0,0,0,0,0;];
E=[ 0,0,0,0,0,0,0,0;
    0,(1-delta)*cap,0,0,0,0,0,0;
    0,alpha,0,0,0,0,0,0;
    0,0,0,0,0,0,0,0;
    0,0,0,0,lambda,0,0,0;
    0,0,0,0,0,1,0,0;
    0,0,0,0,0,0,1,0;
    0,0,0,0,0,0,0,1;];
%use gensys
[G,C,impact,fmat,fwt,ywt,gev,eu]=gensys(D,E,zeros(8,1),[0;0;0;0;1;0;0;0],[0,0,0; 0,0,0; 0,0,0; 0,0,0; 0,0,0; 1,0,0; 0,1,0; 0,0,1]);

F = G([1,2,3,4,5],[1,2,3,4,5]);
%H = [1 0 0 0 0; 0 0 1 0 0; 0 0 0 1 0];
H = [eye(4),zeros(4,1)];%G([1,2,3,4],[2,5]);
Q = vol*impact([1,2,3,4,5])*impact([1,2,3,4,5])';
%R = diag([sigmay,sigmac,sigmal]);
R = diag([sigmay,sigmak,sigmac,sigmal]);
initx = [y0(1);k0(1);c0(1);l0(1);0];
%initx = [y0(1);0;c0(1);l0(1);0];
initV = 0.001*eye(5);

%y1=[y0,c0,l0]';
y1=y2';
[xfilt, Vfilt, VVfilt, loglik] = kalman_filter(y1, F, H, Q, R, initx, initV);

y1=y2';
figure(1)
clf
for i=1:4
subplot(2,2,i)
hold on
plot(time, y1(i,:), 'b*-');
errorbar(time, xfilt(i,:), sqrt(reshape(Vfilt(i, i, :),1,size(Vfilt,3))), 'rx:');
hold off
legend('observed', 'filtered', 2)
xlabel('time')
if i==1
ylabel('output')
elseif i==2
    ylabel('capital')
elseif i==3
    ylabel('consumption')
elseif i==4
    ylabel('labor')
end

end
%[beta, g, ro, alpha, fi ,lambda, delta, sigmay, sigmak, sigmac, sigmal]


figure(2)
clf
hold on
plot(time, y1(1,:), 'b*-');
errorbar(time, xfilt(5,:), sqrt(reshape(Vfilt(5, 5, :),1,size(Vfilt,3))), 'rx:');
hold off
legend('output', 'productivity filtered', 2)
xlabel('time')
ylabel('productivity vs output')

[beta, g, ro, alpha, fi ,lambda, delta, sigmay, sigmak, sigmac, sigmal, vol]


