clear
%!del kprdea.res
%diary kprdea.res
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           Canonical R.B.C. Model            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%**********************
% Calibration and steady state
%**********************



% Technology stochastic process 

rhoa=.95;
sdepsilona=.0072;
%sdepsilona=.0012;
vara=(sdepsilona^2)/(1-rhoa^2);
ETTCHNO=vara^.5;


%production function Y=AK^(alpha)(XN)^(1-alpha)

alpha=.36;
gamma=1.009;


%preferences log C + theta/1-eta L^(1-eta)

N=.2;
eta=0;

%steady state ratio
isk=.025;
ysk=0.12;
ysc=(ysk/(ysk-isk));


delta=isk-(gamma-1);
beta=gamma/(alpha*ysk+1-delta);
theta=(1-alpha)*ysc*((1-N)^eta)/N;


Nsk=ysk^(1/(1-alpha));
k=N*(1/Nsk);
y=k*ysk;
I=k*isk;
c=y-I;
lambda=1/c;
w=(1-alpha)*y/N;
R=gamma/beta+delta-1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Displaying steady state values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

param1=[alpha    beta    delta     N     rhoa     eta];
param2=[y     c     I     k     w    R];
disp('   alpha       beta      delta        N     rhoa    eta'); 
disp(' ');
disp(param1);
disp(' ');
disp('    y          c          I          k           w        R');
disp(param2);

pause;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Linearized Model in Matrix Form
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



M1=zeros(6,6);

M2=zeros(6,3);
M3I=zeros(3,3);
M3L=zeros(3,3);

M4I=zeros(3,6);
M4L=zeros(3,6);

M5=zeros(3,1);

%%%%consumption (1) %%%%

M1(1,1)=-1;
M2(1,3)=1;


%%%%labor supply (2)%%%

M1(2,2)=eta*N/(1-N);
M1(2,5)=-1;
M2(2,3)=1;


%%%%%good market equilibrium (3)%%%%%

M1(3,3)=y;
M1(3,1)=-c;
M1(3,4)=-I;

%%%%%%Production function (4)%%%%%

M1(4,3)=1;
M1(4,2)=-(1-alpha);
M2(4,1)=alpha;
M2(4,2)=1;

%%%Capital demand (5)%%%%

M1(5,3)=1;
M1(5,6)=-1;
M2(5,1)=1;

%%%%%Labor demand (6)%%%

M1(6,3)=1;
M1(6,5)=-1;
M1(6,2)=-1;

%%%%%Keynes-Ramsey (I)%%%%

M3I(1,3)=gamma;
M3L(1,3)=-gamma;
M4I(1,6)=-beta*R;

%%%%% Capital dynamics(II) %%%%

M3I(2,1)=gamma;
M3L(2,1)=-1+delta;
M4L(2,4)=isk;

%%%%Technology stochastic process (III)%%%%

M3I(3,2)=1;
M3L(3,2)=-rhoa;
M5(3,1)=1;

%%%%%%%%%%%%%%
%Reduced Form
%%%%%%%%%%%%%%


L0=M3I-M4I*inv(M1)*M2;
L1=M4L*inv(M1)*M2-M3L;
L2=M5;

W=L0\L1;
Q=L0\L2;

%%%%%%%%%%%%%%%%
%Eigenvalues
%%%%%%%%%%%%%%%%

[PP,MU]=eig(W);

MU=abs(MU);
[llambda,kk]=sort(diag(MU));

P1=PP(:,kk);
P1I=inv(P1);
MU1=P1I*W*P1;



[abf abf] = size(MU1);

ab=0;
for i = 1:abf;
if abs(MU1(i,i)) < 1,
ab=ab+1;
else;
end;
end;

af = abf- ab;

for i = 1:abf;
if abs(MU1(i,i)) == 1,
disp('Unit root')
else;
end;
end;

disp('backward    ')
disp(ab)
disp('forward')
disp(af)

disp('Eigenvalues')
disp(' ')
disp(diag(MU1));

pause



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Saddle Path Condition 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 
%coordonnées /variable backward
P1IB=[P1I(3,1:2)];

%coordonnées / variable forward
P1IF=[-P1I(3,3)];

%expression forward en fonction backward : condition initiale
GG=inv(P1IF)*P1IB;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Policy rules 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%dynamique variable backward
WB=[W(1,1:2)]
WF=[W(1,3)];
PIBE=WB+(WF*GG);
PIBC=[W(2,1:2)];
PIB=[PIBE;PIBC];


%dynamique variable de contrôle
M2B=[M2(1:6,1:2)];
M2F=[M2(1:6,3)];
PIC=M1\(M2B+M2F*GG);


%%%%%%%%%%%%%%%%%%%
%Regle de decisions
%%%%%%%%%%%%%%%%%%%


disp('State Dynamics (k,A)') 
disp(' ')
disp(PIB)


disp('Policy rules (c,H,y,I,w,R) aux etats')
disp(' ')
disp(PIC)


pause


%%%%%%%%%%%%
%Impulse Response Function
%%%%%%%%%%%%

nrep=100;

CHOC=[0 ; 1]

for j=1:nrep;
RC=(PIB^(j-1))*CHOC;
RCK(j)=RC(1);
RCA(j)=RC(2);
end;

for j=1:nrep;
RC=(PIC*PIB^(j-1))*CHOC;
 RCC(j)=RC(1);
 RCH(j)=RC(2);
 RCY(j)=RC(3);
 RCI(j)=RC(4);
 RCW(j)=RC(5);
 RCR(j)=RC(6);

end;

figure
subplot(221),plot(RCK(1:100))
title('Capital')
xlabel('quarters')
ylabel('% Dev.   ')

subplot(222),plot(RCA(1:100))
title('Technology')
xlabel('quarters')
ylabel('% Dev.   ')

subplot(223),plot(RCW(1:100))
title('Labor Productivity')
xlabel('quarters')
ylabel('% Dev.   ')

subplot(224),plot(RCR(1:100))
title(' Capital rental rate ')
xlabel('quarters')
ylabel('% Dev.   ')

pause;

figure
subplot(221),plot(RCC(1:100))
title('Consumption')
xlabel('quarters')
ylabel('% Dev.   ')

subplot(222),plot(RCY(1:100))
title('Output')
xlabel('quarters')
ylabel('% Dev.   ')


subplot(223),plot(RCH(1:100))
title('Hours')
xlabel('quarters')
ylabel('% Dev.   ')

subplot(224),plot(RCI(1:100))
title(' investment')
xlabel('quarters')
ylabel('% Dev.   ')

pause;

%*********************************************
% Stochastic simulation
%*********************************************

nsimul=100;

nlong=160;


for j=1:nsimul;

disp('simulation')
disp(j)

% simulation des jiemes parties cycliques 

aleaa(1:nlong,j)=randn(nlong,1);

% tirage des innovations

for i=1:nlong;
epsa(i)= aleaa(i,j) * (sdepsilona);
end;


% construction des chocs CHT

CHT(1)=epsa(1);

for i=2:nlong;
CHT(i)=rhoa*CHT(i-1)+epsa(i);
end;

% initialisation de la partie cyclique du capital et de l'emploi

KC(1)=0;


% parties cycliques

for i=2:nlong;
KC(i)=PIB(1,1)*KC(i-1)+PIB(1,2)*CHT(i-1);

end;

for i=1:nlong;

YC(i)=PIC(3,1)*KC(i)+PIC(3,2)*CHT(i);

CC(i)=PIC(1,1)*KC(i)+PIC(1,2)*CHT(i);

NC(i)=PIC(2,1)*KC(i)+PIC(2,2)*CHT(i);

IC(i)=PIC(4,1)*KC(i)+PIC(4,2)*CHT(i);

PC(i)=PIC(5,1)*KC(i)+PIC(5,2)*CHT(i);

end;



invyc=YC*YC';

pmc(j)=invyc\(YC*CC');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Business Cycles Properties
%***********

% HP FILTER matrix

% lambda of the Hp filter
ntrunc=1;
lhp=1600;

MHP=[1+lhp -2*lhp lhp zeros(1,nlong-ntrunc+1-3);...
-2*lhp 1+5*lhp -4*lhp lhp zeros(1,nlong-ntrunc+1-4);...
zeros(nlong-ntrunc+1-4,nlong-ntrunc+1);...
zeros(1,nlong-ntrunc+1-4) lhp -4*lhp 1+5*lhp -2*lhp;...
zeros(1,nlong-ntrunc+1-3) lhp -2*lhp 1+lhp];

for i=3:nlong-ntrunc+1-2
MHP(i,i-2)=lhp;
MHP(i,i-1)=-4*lhp;
MHP(i,i)=1+6*lhp;
MHP(i,i+1)=-4*lhp;
MHP(i,i+2)=lhp;
end

%Filtering the series

YHP=YC'-MHP\YC';
NHP=NC'-MHP\NC';
CHP=CC'-MHP\CC';
IHP=IC'-MHP\IC';
PHP=PC'-MHP\PC';



ETNHP(j)=std(NHP(1:nlong));
ETYHP(j)=std(YHP(1:nlong));
ETCHP(j)=std(CHP(1:nlong));
ETIHP(j)=std(IHP(1:nlong));
ETPHP(j)=std(PHP(1:nlong));

RHOCHP(j)=CORR(YHP,CHP);
RHONHP(j)=CORR(YHP,NHP);
RHOIHP(j)=CORR(YHP,IHP);
RHOPHP(j)=CORR(YHP,PHP);

RHOPNHP(j)=CORR(NHP,PHP);


xxN=NHP(2:nlong);
yyN=NHP(1:nlong-1);
RHONNHP(j)=CORR(xxN,yyN);

xxC=CHP(2:nlong);
yyC=CHP(1:nlong-1);
RHOCCHP(j)=CORR(xxC,yyC);

xxY=YHP(2:nlong);
yyY=YHP(1:nlong-1);
RHOYYHP(j)=CORR(xxY,yyY);

xxI=IHP(2:nlong);
yyI=IHP(1:nlong-1);
RHOIIHP(j)=CORR(xxI,yyI);

xxP=YHP(2:nlong);
yyP=YHP(1:nlong-1);
RHOPPHP(j)=CORR(xxP,yyP);



end;



mETNHP=mean(ETNHP);
mETYHP=mean(ETYHP);
mETCHP=mean(ETCHP);
mETIHP=mean(ETIHP);
mETPHP=mean(ETPHP);

mRHOCHP=mean(RHOCHP);
mRHONHP=mean(RHONHP);
mRHOIHP=mean(RHOIHP);
mRHOPHP=mean(RHOPHP);
mRHOPNHP=mean(RHOPNHP);

mRHOCCHP=mean(RHOCCHP);
mRHOYYHP=mean(RHOYYHP);
mRHONNHP=mean(RHONNHP);
mRHOIIHP=mean(RHOIIHP);
mRHOPPHP=mean(RHOPPHP);


mET1=[mETNHP mETYHP mETCHP mETIHP mETPHP];
mET1R=[mETNHP mETYHP mETCHP mETIHP mETPHP]./mETYHP;
mRHO1=[mRHONHP 1 mRHOCHP mRHOIHP mRHOPHP];
mARHO1=[mRHONNHP mRHOYYHP mRHOCCHP mRHOIIHP mRHOPPHP];

disp('variable order: N - Y - C - I - P')
disp(' ')

disp('standard deviation')
disp(' ')
disp(mET1)



disp('relative standard deviation')
disp(' ')
disp(mET1R)


disp('correlation')
disp(' ')
disp(mRHO1)


disp('serial correlation')
disp(' ')
disp(mARHO1)

disp('correlation N-W')
disp(' ')
disp(mRHOPNHP)


break

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Calcul du welfare

WW=log(c)+eta*log(1-H)-1/(2*(c)^2)*( c^2*(mETCHP)^2 )-1/2*eta/(1-H)^2*H^2*(mETNHP)^2;

WW1= log(c) + eta* log(1-H) ;

WW2=-1/(2*(c)^2)*( c^2*(mETCHP)^2 )-1/2*eta/(1-H)^2*H^2*(mETNHP)^2;



% expression du welfare total en point de conso par rapport au %steady-state de r‚f‚rence


CC=100 * (( 1/c * ( exp(WW - eta*log(1-H)) )) - 1 );



disp('     W      W steady  W variance Pts de conso')
disp([WW  WW1 WW2 CC])