%RAMSEY_MOMS.M
%table = ramsey_moms(model) computes a table of second moments implied by the Ramsey equilibrium 
%in the model developed by ``Comparing Two Variants of Calvo-Type Wage Stickiness,'' by 
%Stephanie Schmitt-Grohe and Martin Uribe (2006).
%Layout of the Table:
%Column 1: Standard Deviations
%Column 2: Serial Correlations
%Column 3: Correlation with Output Growth
%Row 1: Output Growth 
%Row 2: Consumption Growth 
%Row 3: Investment Growth 
%Row 4: Nominal Interest Rate
%Row 5: Price Inflation
%Row 6: Wage Inflation
%Row 7: Hours
%Input: 
%model = a string variable taking the values 
%'ehl' or 'sgu' indicating the variant of wage inertia to be used. 
%The default value is 'sgu'.
%(c) Stephanie Schmitt-Grohe and Martin Uribe, October 2006.




function table = ramsey_moms(model)

if nargin<1
    model = sgu
end

%Compute first-order approximation to the policy funciton of the Ramsey problem
[gx, hx, controlvar_cu, statevar_cu, eta] = ramsey_run(model);


nx=size(hx,1);
ny=size(gx,1);
nv = nx+ny; %total number of variables (states+controls) in the system

%Symbolic vector of control and state varaibles
var_cu = [controlvar_cu statevar_cu];

%Index of variables in vector var_cu
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');
niv = find(var_cu=='iv_cu');
nw = find(var_cu=='w_cu');
ns = find(var_cu=='s_cu');
nstil = find(var_cu=='stil_cu');
nk = find(var_cu=='k_cu');
ng = find(var_cu=='g_cu');
nmuz = find(var_cu=='muz_cu');
nmuupsilon = find(var_cu=='muupsilon_cu');
nu = find(var_cu=='u_cu');
nvt = find(var_cu=='vt_cu');


%Matrix Scaling Standard Deviations of Innovation to State Vector
varshock = eta*eta';

%Let the vector z_t = [y_t;x_t] Then, up to 1st order 
%z_t = Az x_t + Bz eps_{t+1}
Az = [gx;eye(nx)];

Bz = zeros(nv,size(eta,2));

%z_{t+1} =Azp x_t+Bzp epx_{t+1}
Azp = Az * hx;
Bzp = Az * eta;

%output
Aoutput = Az(noutput,:);
Boutput = Bz(noutput,:);


%Hours
Ahdp = Azp(nhd,:);
Bhdp = Bzp(nhd,:);

%outputp
Aoutputp = Azp(noutput,:);
Boutputp = Bzp(noutput,:);


%Consumption
Ac = Az(nc,:);
Bc = Bz(nc,:);

%Future consumption
Acp = Azp(nc,:);
Bcp = Bzp(nc,:);

%Investment
Aiv = Az(niv,:);
Biv = Bz(niv,:);

%Future investment
Aivp = Azp(niv,:);
Bivp = Bzp(niv,:);


%Real wage
Aw = Az(nw,:);
Bw = Bz(nw,:);

%Future real wage
Awp = Azp(nw,:);
Bwp = Bzp(nw,:);


%Nominal Interest Rate
Arp = Azp(nr,:);
Brp = Bzp(nr,:);


%Price Inflation
Apai = Az(npai,:);
Bpai = Bz(npai,:);


%Future Price Inflation
Apaip = Azp(npai,:);
Bpaip = Bzp(npai,:);


%Investment-Specific Tech Shock
Amuupsilonp = Azp(nmuupsilon,:);
Bmuupsilonp = Bzp(nmuupsilon,:);

%Neutral Tech. Shock
Amuzp = Azp(nmuz,:);
Bmuzp = Bzp(nmuz,:);

%load parameter THETA
eval(['load ramsey_ss_' model ' THETA'])

%muzstar = (THETA/(1-THETA)) muupsilon + muz;
Amuzstarp = (THETA/(1-THETA)) * Amuupsilonp + Amuzp;
Bmuzstarp = (THETA/(1-THETA)) * Bmuupsilonp + Bmuzp;

%Output Growth
Agoutputp = Aoutputp - Aoutput + Amuzstarp;
Bgoutputp = Boutputp - Boutput + Bmuzstarp;

%Consumption Growth
Agcp = Acp - Ac + Amuzstarp;
Bgcp = Bcp - Bc + Bmuzstarp;


%Investment Growth
Agivp = Aivp - Aiv + Amuzstarp + Amuupsilonp;
Bgivp = Bivp - Biv + Bmuzstarp + Bmuupsilonp;

%Real Wage Growth
Agwp = Awp - Aw + Amuzstarp;
Bgwp = Bwp - Bw + Bmuzstarp;

%Wage Inflation
Apaiwp = Agwp + Amuzstarp + Apaip;
Bpaiwp = Bgwp + Bmuzstarp + Bpaip;

%xi_{t+1} = [goutputp;gcp;givp;hdp;rp;paip;paiwp]; Then xi_{t+1}=A*x_t + B*eps_{t+1}, where
A = [Agoutputp;Agcp;Agivp;Arp;Apaip;Apaiwp;Ahdp];

B = [Bgoutputp;Bgcp;Bgivp;Brp;Bpaip;Bpaiwp;Bhdp];

%Compute var-cov matrix of x_t
[VARY,VARX]=mom(gx,hx,varshock,0);

%Compute Var-Cov matrix of XI_t
VARXI = A*VARX*A'+B*B';

table=zeros(7,3);
table(:,1)=400*sqrt(diag(VARXI));
table(7,1) = 100*sqrt(VARXI(end,end));

table(:,3) = VARXI(:,1)./sqrt(diag(VARXI))/sqrt(VARXI(1,1));


%Compute 1st order -Cov matrix of x_t with x_{t-1}
[VARY,VARX]=mom(gx,hx,varshock,1);

%Compute 1st order  Cov matrix of XI_t with XI_{t-1}
VAR1XI = A*VARX*A'+ A*eta*B';
table(:,2)=(diag(VAR1XI))./diag(VARXI)

disp('Columns: s.d, serial correlation, correlation with output growth.')
disp('Rows: output growth, consumption growth, investment growth, nominal interest rate, price inflation, wage inflaiton, hours.')