%Econ 202C Macro. Homework#2. Problem1 Anton Cheremukhin
close all;  clear all;  clc;
%data
[time, gdp, cons, inv1, hours,inv2] = textread('datahw2.txt','%d %f %f %f %f %f');
figure(1);
plot(time,gdp,time,cons,time,inv1,time,hours)
legend('gdp','cons','inv','hours')
figure(2);
plot(time,cons./gdp,time,inv1./gdp,time,inv2./gdp)
legend('c/y','inv/y','(inv+dur)/y')
w=((cons./gdp).*hours)/((cons(1)/gdp(1))*hours(1));
T=size(w,1);
tau=log(w);
figure(3);
plot(time,tau)
legend('wedge')
Y=tau(2:T);  X=tau(1:T-1);
ro=inv(X'*X)*(X'*Y);
v=Y-X*ro;

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;
i0=log(inv1);   X=[ones(T,1) time]; beta=inv(X'*X)*(X'*i0); i0=i0-X*beta;
l0=log(hours);   X=[ones(T,1) time]; beta=inv(X'*X)*(X'*l0); l0=l0-X*beta;

figure(4);
plot(time,y0,time,c0,time,i0,time,l0)
legend('gdp','cons','inv','hours');

D=[ 0,-1,-1/3,1,0,0;
    8/53,-1,0,0,1,-8/53;
    50/36,3/4,-2/3,0,0,0;
    0,0,0,1,0,0;
    0,1,0,0,0,0;
    0,0,1,0,0,0];
E=[ -1/3,0,0,0,0,0;
    0,0,0,0,0,0;
    41/36,0,0,0,0,0;
    0,0,0,ro,0,0;
    0,0,0,0,1,0;
    0,0,0,0,0,1];
%use gensys
[G,C2,impact2,fmat,fwt,ywt,gev,eu]=gensys(D,E,zeros(6,1),[0;0;0;1;0;0],[0,0;0,0;0,0;0,0;1,0;0,1])
%technology shocks definition
N=T-1;
%generate sequence for t=1:N
%z=sigma*randn(1,N);
z=v;
%z=zeros(1,N);  %z(10)=2*sigma;
G2=G(1:4,1:4);
x=zeros(4,N);
for i=2:N
    x(:,i)=G2*x(:,(i-1))+impact2(1:4)*z(i);
end
yt=(x(1,:)+2*x(3,:))/3;
it=(50/36)*x(1,2:T-1)-(41/36)*x(1,1:T-2);

%simulated data
A=corrcoef([x(2:3,1:45)' yt(1:45)' it'])
%real data
B=corrcoef([c0 l0 y0 i0])

%plot results
t=1:1:N;
figure(5)
subplot(2,2,1);
plot(t,y0(2:N+1),'b',t,yt(1:N),'r')
legend('output','output sim')
subplot(2,2,2);
plot(t,c0(2:N+1),'b',t,x(2,1:N),'r')
legend('consumption','consumption sim')
subplot(2,2,4);
plot(t,l0(2:N+1),'b',t,x(3,1:N),'r')
legend('labor','labor sim')
t=1:1:N-1;
subplot(2,2,3);
plot(t,i0(2:N),'b',t,it(1:N-1),'r')
legend('investment','investment sim')
t=1:1:N;
figure(6)
plot(t,w(2:N+1)-1,'b',t,x(4,1:N),'r')
