%RAMSEY_MOMENTS.M
%Produce a table contaning first and second moments of selected variables in the Ramsey equilibrium under income taxation for the model in ``Optimal Fiscal and Monetary Policy in a Medium-Scale Macroeconomic Model,'' by Stephanie Schmitt-Grohe and Martin Uribe, April 2005 
%(c) April 5, 2005 by Stephanie Schmitt-Grohe and Martin Uribe

%Compute analytically the Ramsey equilibrium conditions and its first derivatives, evaluate this functions at the Ramsey steady state, and compute matrices gx and hx defining a first-order approximation to the Ramsey equilbrium

ramsey_dynamics

disp('execution of RAMSEY_DYNAMICS.M completed')


%Number of states
nx=size(hx,1);

%Number of controls
ny=size(gx,1);

%Number of variables
nv=nx+ny;

T=200; %Simulation horizon

n = 1000; %number of simulated economies

%VarCov matrix of exog. shocls
std_matrix = diag([STD_EPSN; STD_EPSG; STD_EPSZ]);

%Set initial seed
randn('state',sum(100*clock));

%all variables 
var_cu = [controlvar_cu statevar_cu];

%Index of variables
nqq = find(var_cu=='qq_cu');
nhd = find(var_cu=='hd_cu');
noutput = find(var_cu=='output_cu');
nc = find(var_cu=='c_cu');
nla = find(var_cu=='la_cu');
npai = find(var_cu=='pai_cu');
nrk = find(var_cu=='rk_cu');
nf2 = find(var_cu=='f2_cu');
nx2 = find(var_cu=='x2_cu');
nr = find(var_cu=='r_cu');
ntauh = find(var_cu=='tauh_cu');
ntauk = find(var_cu=='tauk_cu');
niv = find(var_cu=='iv_cu');
nw = find(var_cu=='w_cu');
ns = find(var_cu=='s_cu');
nstil = find(var_cu=='stil_cu');
nassets = find(var_cu=='assets_cu');
nk = find(var_cu=='k_cu');
nn = find(var_cu=='n_cu');
ng = find(var_cu=='g_cu');
nz = find(var_cu=='z_cu');

npick = [ntauk nr npai noutput nc niv nhd nw nassets];

%initialize table
table = zeros(numel(npick),4);

for jn=1:n

%Initial condition
x0 = zeros(nx,1);

%Generate innovations
shocks = zeros(T,nx);
shocks(:,end-2:end) = sign(randn(200,3)) * std_matrix;

%Generate path of endog. variables

for j=1:T

%Path of controls
Y(j,1:ny) = (gx*x0)';

%Path of states
X(j,1:nx) = x0';

x0 = hx*x0 + shocks(j,:)';

end %for j=1:T

V = [Y X];

stdV = std(V);

corrV = corrcoef(V);

corrVoutput = corrV(:,noutput);

corrcoef([V(2:end,:) V(1:end-1,:)]);

sercorrV = diag(ans(1:nv,nv+1:end));

nvar = [ncontrolvar;nstatevar];

table_jn = [nvar(npick) stdV(npick)' sercorrV(npick) corrVoutput(npick)]; 

weight_jn = 1/jn;

table = table * (1-weight_jn) + table_jn * weight_jn;

end  %for jn

table(1,1) = table(1,1)*100;

table(2,1)=(table(2,1)^4-1)*100;

table(3,1)=(table(3,1)^4-1)*100;

table(:,2)=table(:,2)*100;

table(2:3,2) = table(2:3,2)*4;

table(:,1:4) = round(table(:,1:4)*100)/100;

table

disp('NOTE')

disp('Columns: (1) mean, (2) std. dev., (3) ser. corr., (4) corr w. output')

disp('Rows: (1) tax rate, (2) nominal interest rate, (3) inflation, (4) output, (5) consumption, (6) investment, (7) hours, (8) real wage, (9) real government liabilities')