?? rbc_eng.m
字號:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% The baseline RBC model %
% %
% %
% U(C,1-h)=log(C)+theta log(1-h) %
% %
% Y=a(K^alpha)(h^(1-a)) %
% %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; % Clear the memory
delete rbc.res; % If it exists, delete the file exo.res
diary rbc.res; % Open a new file called exo.res that will store
% all what is written on the screen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Algorithm Parameters %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ncont=5; % # of static equations
nback=1; % # of Backward endogenous state variables
nshoc=1; % # of exogenous shocks
ntotb=nback+nshoc; % Total # of Backward state variables
nforw=1; % # of Forward endogenous state variables
nstat=ntotb+nforw; % Total # of state variables
long=120; % Length of simulated series
trunc=50; % Truncature of simulations
slong=long+trunc; % Length of simulation
nsim=5000; % # of simulations
nrep=20; % IRF horizon
Mcc=zeros(ncont,ncont); % Matrix Control-Control
Mcs=zeros(ncont,nstat); % Matrix Control-State
Mss0=zeros(nstat,nstat); % Matrix State-State
Mss1=zeros(nstat,nstat); % Matrix State-State lagged
Msc0=zeros(nstat,ncont); % Matrix State-Control
Msc1=zeros(nstat,ncont); % Matrix State-Controle lagged
Mse=zeros(nstat,nshoc); % Matrix State-Shocks
select=[1:5]; % Selection of variables of interest
indy=3; % Index of Y(t) in the selection
indn=2; % Index of h(t) in the selection
indp=5; % Index of Y/N(t) in the selection
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Structural Parameters of the Economy %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The model is calibrated for the US econmy on Quarterly data
% The numbers are taken from Cooley and Prescott (1995) Chapter 1
%
whoy=0.6; % Share of wages in total value added
ykoy=3.32; % Yearly K/Y
yiok=0.076; % Yearly I/K
gzy=0.0156; % Yearly per capita technological rate of growth
gny=0.0120; % Yearly population rate of growth
hss=0.31; % People devote 31% of their total time endowment to work
rhoa=0.95; % Persistence of technological shock
stda=0.0079; % Standard deviation of technological shock
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Steady state %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gamy=(1+gzy)*(1+gny); % Quarterly total rate of growth
gam=gamy^(1/4); % Quarterly total rate of growth
alpha=1-whoy; % elasticity of output w.r. to capital
dy=1+yiok-(1+gzy)*(1+gny); % Yearly depreciation rate
by=gamy/(alpha*(1/ykoy)+1-dy); % Yearly discount factor
delta=1-(1-dy)^(1/4); % Quarterly depreciation rate
beta=by^(1/4); % Quarterly discount factor
yok=(gam-beta*(1-delta))/(alpha*beta); % Model based Y/K (quarterly)
iok=gam+delta-1; % Model based I/K (quarterly)
ioy=iok/yok; % Model based I/Y (quarterly)
coy=1-ioy; % Model based C/Y (quarterly)
theta=(1-hss)*(1-alpha)/(hss*coy); % Parameter of the disutility of leisure
yss=(yok^(alpha/(alpha-1)))*hss; % Level of output
css=coy*yss; % Level of consumption
iss=ioy*yss; % Level of investment
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Matrices Coefficients %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Static equations (c h y i p)
%
%
% Consumption
%
cc=-1;
cl=1;
%
% Hours
%
hh=1/(1-hss);
hy=-1;
hc=1;
%
% Production function
%
yy=1;
yh=alpha-1;
yk=alpha;
ya=1;
%
% Investment
%
ic=-coy;
iy=1;
ii=-ioy;
%
% Productivity
%
py=-1;
ph=1;
pp=1;
%
% Dynamic equations (k a lambda)
%
%
% Accumulation
%
kk=gam;
kkl=delta-1;
kil=gam+delta-1;
%
% Shock
%
aa=1;
aal=-rhoa;
ae=1;
%
% Euler
%
ll=-1;
lll=1;
lk=alpha*beta*yok/gam;
ly=lk;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Define the matrices %
% %
% %
% Mcc X(t) = Mcs S(t) %
% Mss0 S(t+1) + Mss1 S(t) = Msc0 X(t+1) + Msc1 X(t) + e(t+1) %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Static equations :
%
% c h y i p
% 1 2 3 4 5
%
% Dynamic equations :
%
% k a lambda
% 1 2 3
%
% Mcc
%
Mcc(1,1)=cc;
Mcc(2,1)=hc;
Mcc(2,2)=hh;
Mcc(2,3)=hy;
Mcc(3,2)=yh;
Mcc(3,3)=yy;
Mcc(4,1)=ic;
Mcc(4,3)=iy;
Mcc(4,4)=ii;
Mcc(5,2)=ph;
Mcc(5,3)=py;
Mcc(5,5)=pp;
%
% Mcs
%
Mcs(1,3)=cl;
Mcs(3,1)=yk;
Mcs(3,2)=ya;
%
% Mss0
%
Mss0(1,1)=kk;
Mss0(2,2)=aa;
Mss0(3,1)=lk;
Mss0(3,3)=ll;
%
% Mss1
%
Mss1(1,1)=kkl;
Mss1(2,2)=aal;
Mss1(3,3)=lll;
%
% Msc0
%
Msc0(3,3)=ly;
%
% Msc1
%
Msc1(1,4)=kil;
%
% Mse
%
Mse(2,1)=ae;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Solve the Farmer's system %
% %
% -1 %
% X(t) = Mcc Mcs S(t) %
% %
% -1 -1 %
% S(t+1)=(Mss0-Msc0 Mcc Mcs)(Msc1 Mcc Mcs-Mss1)S(t) %
% -1 %
% +(Mss0-Msc0 Mcc Mcs)e(t+1) %
% %
% S(t+1) = W S(t) + R e(t+1) %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% First form the preceeding system
%
M0=inv(Mss0-Msc0*inv(Mcc)*Mcs);
M1=(Mss1-Msc1*inv(Mcc)*Mcs);
W=-M0*M1;
%
% 1) Compute the eigenvalues (MU) and eigenvectors (P) of W
% 2) Compute their modulus (AMU)
% 3) Sort the eigenvalues and rearrange the matrix of eigenvectors accordingly (P)
% 4) Form the inverse of P (Ps)
%
[P,MU] = eig(W);
AMU=diag(abs(MU))';
disp('Eigenvalues / modulus of the system :')
disp('=====================================')
disp(' ');
disp([diag(MU)';AMU])
flag=sum(AMU>1);
if flag==nforw;
disp('Blanchard-Kahn conditions checked')
else
disp('Blanchard-Kahn conditions NOT checked')
end;
disp(' ');
[MU,k] = sort(AMU);
P=P(:,k);
Ps=inv(P);
%
% Solves the system using Farmer's method
%
M=[eye(ntotb);-inv(Ps(ntotb+1:nstat,ntotb+1:nstat))*Ps(ntotb+1:nstat,1:ntotb)];
MSS=W(1:ntotb,:)*M;
M2=M0*Mse;
MSE=M2(1:ntotb,nshoc);
PI=inv(Mcc)*Mcs*M;
%
% We are done !!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% So we now get : %
% %
% X(t)=Pi S(t) %
% %
% S(t+1)=MSS S(t) + MSE e(t+1) %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Policy Functions :')
disp('==================');
disp(' ');
disp('K(t+1) - A(t+1) as functions of (K(t),A(t)) :');
disp(' ');
disp(MSS);
disp('C(t) - h(t) - Y(t) - I(t) - Y(t)/h(t) as functions of (K(t),A(t)) :');
disp(' ');
disp(PI);
diary off; % Stop recording output in a file
disp('Hit a Key');pause;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Impulse response functions to a technological shock %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CA=[1]; % initial value of the impulse
REPSA=zeros(ntotb,nrep); % IRF matrix of state variables
REPSA(:,1)=MSE*CA; % Initialize the IRF matrix of state variables
for i=2:nrep; %
REPSA(:,i)=MSS*REPSA(:,i-1); % Main loop for IRF
end; %
REPXA=PI*REPSA; % IRF matrix of controls
%
% Graphs
%
T=1:nrep;
subplot(221);plot(T,REPXA(1,:));title('IRF(C,A)');
xlabel('Quarters');ylabel('% deviation')
subplot(222);plot(T,REPXA(3,:));title('IRF(Y,A)');
xlabel('Quarters');ylabel('% deviation')
subplot(223);plot(T,REPXA(4,:));title('IRF(I,A)');
xlabel('Quarters');ylabel('% deviation')
subplot(224);plot(T,REPSA(1,:));title('IRF(K,A)');
xlabel('Quarters');ylabel('% deviation')
print -dps irbc1.eps;
pause;close % wait for you to hit a key and close the graph window
subplot(221);plot(T,REPXA(2,:));title('IRF(H,A)');
xlabel('Quarters');ylabel('% deviation')
subplot(222);plot(T,REPXA(5,:));title('IRF(P,A)');
xlabel('Quarters');ylabel('% deviation')
print -dps irbc2.eps;
pause;close % wait for you to hit a key and close the graph window
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Simulation %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
HP=hpmat(long,1600); % Retrieve the matrix for the HP filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Loop of Simulation %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The deterministic component
%
trend=(1:long)*log(gam);
determ=[1;0;1;1;1]*trend+log([css;hss;yss;iss;yss/hss]*ones(size(trend)));
%
% stochastic Simulation
%
SSIM=zeros(ntotb,slong); % matrix of simulated state variables
SHP=zeros(nsim,length(select)); % matrix to store simulated standard deviations
SRHP=zeros(nsim,length(select)); % matrix to store simulated relative standard deviations
CHP=zeros(nsim,length(select)); % matrix to store simulated correlations
CHYN=zeros(nsim,1); % matrix to store simulated employment-productivity correlation
%
% Beginning of Monte-Carlo loop
%
for j=1:nsim;
disp(sprintf('Simulation : %4.0f',j));
SHOCK=[stda]*randn(nshoc,slong); % Draw a random number from N(0,stda)
SSIM=MSE*SHOCK(1,1); % Initialize each simulation
for t=2:slong; %
SSIM(:,t)=MSS*SSIM(:,t-1)+MSE*SHOCK(:,t); % Simulate state variables
end; %
%
% Retrieve control variables. Note that :
% 1) we add the deterministic component
% 2) we throw out the first trunc simulations in order
% to be free of any initial conditions effects
%
XSIM=PI(select,:)*SSIM(:,trunc:trunc+long-1)+determ;
%
% HP-filter the control variables
%
HPSIM=XSIM'-HP\XSIM';
%
% Compute and store the moments of interest
%
SHP(j,:)=std(HPSIM);
SRHP(j,:)=SHP(j,:)/std(HPSIM(:,indy));
VCOV=cov(HPSIM);
CHP(j,:)=VCOV(indy,:)./sqrt(diag(VCOV)'*VCOV(indy,indy));
CHYN(j,:)=VCOV(indn,indp)./sqrt(VCOV(indn,indn)*VCOV(indp,indp));
end;
%
% End of Monte-Carlo loop
%
%
% Now, display the results
%
diary rbc.res; % restart recording output in the file rbc.res
disp(' ')
disp('Quantitative Evaluation :');
disp('=========================');
disp(' ');
disp(sprintf('Based on %5.0f simulations',nsim));
disp(' ');
disp('Standard Deviation : C - h - Y - I - Y/N');
disp(' ');
disp(100*[mean(SHP);std(SHP)]);
disp('Relative Standard Deviation : C - h - Y - I - Y/N');
disp(' ');
disp([mean(SRHP);std(SRHP)]);
disp('Correlation with Output : C - h - Y - I - Y/N');
disp(' ');
disp([mean(CHP);std(CHP)]);
disp('Correlation Hours - productivity');
disp(' ');
disp([mean(CHYN);std(CHYN)]);
diary off;
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -