// DYNARE illustration of 
// T. Sargent, N. Williams and T. Zha (2005) "Shocks and Government Beliefs: The Rise and Fall of American Inflation".
// The paramter values are obtained from simulate_swz.m by Tao Zha, September 2005.
// Michel Juillard, October 2005.
//

//
// suppress warnings
//
warning off;

// 
// endogenous variables:
// pi: inflation rate
// u:  unemployment rate
// x:  government instrument (one-period ahead conditional expectation of inflation)
// L1: Lagrange multiplier on (pi_t = x_{t-1}+sigma2*w2)
// L2: Lagrange multiplier on (u_t = a0*pi_t+a1*pi_{t-1}+a2*u_{t-1}+a3*pi_{t-2}+a4*u_{t-2}+a5+sigma*w)
// 
var L1 L2 pi u x ;
//
// exogenous variables:
// w:  shock on misspecified Phillips curve
// w1: shock on true Phillips curve
// w2: shock on inflation
//
varexo w w1 w2;

//
// parameters of the government optimization problem based on the misspecified Phillips curve
//
parameters pistar ustarstar theta0 theta1 tau1 sigma sigma2 delta lambda a0 a1 a2 a3 a4 a5;
pistar = 2;
ustarstar = 1;
lambda = 1;
delta = 0.9936;
sigma = sqrt(35.653810749685604);
sigma2 = sqrt(18.976712248466477);

//
// Parameters of the true Phillips curve
//
ustar = 6.1104238781343128e+000;
theta0 = -8.8446077428639834e-004;
theta1 = -1.2185510508522781e-002;
tau1 = 9.8917842808371093e-001;
sigma1 = sqrt(35.653810749685604);

//
// governement optimization problem with misspecified Phillips curve
//
model(linear);
u = a5+a0*pi+a1*pi(-1)+a2*u(-1)+a3*pi(-2)+a4*u(-2)+sigma*w;
// Dynare insists on having all variable present in the model at the current period ...
// therefore the 0*x term in next equation
pi = 0*x+x(-1)+sigma2*w2;
2*(pi-pistar)-L2*a0-delta*L2(+1)*a1-delta^2*L2(+2)*a3 +L1 = 0;
2*lambda*(u-ustarstar)+L2-delta*L2(+1)*a2-delta^2*L2(+2)*a4 = 0;
L1(+1) = 0;
end;

shocks;
var w; stderr 1;
var w1; stderr 1;
var w2; stderr 1;
end;

//
// load data for initialization
//
load datainp_upce.prn;
ut = datainp_upce(:,1);
pcet = datainp_upce(:,2);
//
// Let's free some memory
//
clear datainp_upce;      

//
// setting some parameters
//
// maximum lag length (using internal Dynare variable)
maxnlags = ykmin_;
// Effective sample size (excluding 2 lags).
n = length(ut)-maxnlags;

//
// make space for results
//
utpred = zeros(n,1);
xtm1 = zeros(n,1);

//
// Initial value of the coefficient of the misspecified Phillips curve
//
z10 = [
 -1.3236567947738001e-001  1.4191629076828000e-001  1.0927537093697399e+000  -2.1570770229959998e-002  -1.3378436709391001e-001  2.1897028305873001e-001
   ]';
//
// Initialization of the Kalman filter (used only with learning on)
//
P10 = [
 1.0870488505056409e+001  1.4323621422195412e+001  2.2518190299428884e+000  -2.5403728156140705e+001  -9.2792811323474544e-001  -1.0154752905165562e+001
 1.4323621422195412e+001  1.9372095725285959e+001  2.9623869975360826e+000  -3.3983165867677663e+001  -1.1883450541284408e+000  -1.3592338677442111e+001
 2.2518190299428884e+000  2.9623869975360826e+000  4.6904565623517674e-001  -5.2628895615435400e+000  -1.9281142974224980e-001  -2.1050448465146596e+000
 -2.5403728156140705e+001  -3.3983165867677663e+001  -5.2628895615435400e+000  5.9899679355282281e+001  2.1339360050403391e+000  2.3955102199197551e+001
 -9.2792811323474544e-001  -1.1883450541284408e+000  -1.9281142974224980e-001  2.1339360050403391e+000  8.1608300428139130e-002  8.5260055588578354e-001
 -1.0154752905165562e+001  -1.3592338677442111e+001  -2.1050448465146596e+000  2.3955102199197551e+001  8.5260055588578354e-001  9.5810110741206937e+000
    ];
V = [
 8.2322662111882163e+000  -7.7781439991475825e+000  9.2077629698995767e-001  4.9782253939194820e+000  -8.1360515245574161e-001  -4.1414085039716063e+001
 -7.7781439991475825e+000  8.1400447518276913e+000  3.0371661961513308e-002  -5.0890194213372792e+000  1.9352850151333634e+000  6.8591133058425442e+001
 9.2077629698995767e-001  3.0371661961513308e-002  2.9853986140270035e+000  1.1870299279361121e-001  3.7011949370750830e+000  7.2067403473981301e+001
 4.9782253939194820e+000  -5.0890194213372792e+000  1.1870299279361121e-001  3.2031654646433267e+000  -1.0547748262362802e+000  -3.9963285100898247e+001
 -8.1360515245574161e-001  1.9352850151333634e+000  3.7011949370750830e+000  -1.0547748262362802e+000  5.1362258933213196e+000  1.0063609330243516e+002
 -4.1414085039716063e+001  6.8591133058425442e+001  7.2067403473981301e+001  -3.9963285100898247e+001  1.0063609330243516e+002  2.5883091086216796e+003
      ];

zt_tm1 = z10;
Pt_tm1 = P10;
xtm1(1) = pcet(maxnlags+1);
utpred(1) = ustar + tau1*(ut(maxnlags)-ustar);

//
// Index of variables in steady state vector
//
ip = var_index('pi');
iu = var_index('u');
ix = var_index('x'); 
//
// Sargent and William's normalization
//
sig2 = 1.0/sigma1^2;
//
// starting the clock 
//
tic;
//
// looping over observations
//
for ti=2:n
//
// ti1 is index in data vectors with a maxnlags=2 offset
//
   ti1 = ti+maxnlags;
//  
// variable vector for misspecified Phillips curve
//
   phi = [pcet(ti1-1) pcet(ti1-2) ut(ti1-2) pcet(ti1-3) ut(ti1-3) 1];
//  
// variable vector for true Phillips curve
//
   phi2=[phi(1) ut(ti1-1)-ustarstar phi([2 3 6])];

//
// setting model parameters
//
   a0 = zt_tm1(1);
   a1 = zt_tm1(2);
   a2 = zt_tm1(3);
   a3 = zt_tm1(4);
   a4 = zt_tm1(5);
   a5 = zt_tm1(6);
//
// calling Dynare solution routine
// 
   options_.order = 1;
   options_.dr_algo = 0;
   [dr, info] = resol(zeros(endo_nbr,1),0);

//
// building the state vector (the order of variables is the one used in swz2.mod output)
// must be a column vector
//
   ys = dr.ys;
   svec = [pcet(ti1-2)-ys(ip); ut(ti1-2)-ys(iu); xtm1(ti-1)-ys(ix);...
           pcet(ti1-3)-ys(ip); ut(ti1-3)-ys(iu)];
//
// recovering government instrument (and one-period ahead inflation forecast)
// xtm1(ti) = E_{ti-1} pi_{ti}
//
   ixs = var_state_index('x',dr.order_var);
   ius = var_state_index('u',dr.order_var);
//
// standardized inflation prediction error in ti-1
//
   ep = (pcet(ti1-1)-xtm1(ti-1))/sigma2;
//
// standardized unemployment error in ti-1,
// consistent with the government misspecified model
//
   up = (ut(ti1-1)-a5-a0*pcet(ti1-1)-a1*pcet(ti1-2)-a2*ut(ti1-2)-a3*pcet(ti1-3)-a4*ut(ti1-3))/sigma;
   xtm1(ti) = dr.ys(ix)+dr.ghx(ixs,:)*svec+dr.ghu(ixs,:)*[up; 0; ep];
//
// recovering unemployment one-period ahead unemployment forecast
// utpred(ti) = E_{ti-1} u_{ti}
//
   utpred(ti) = ustar + theta0*(pcet(ti1)-xtm1(ti)) + theta1*(pcet(ti1-1)-xtm1(ti-1))...
                + tau1*(ut(ti1-1)-ustar);

//
// Kalman filter updating: replaces zt_tm1 and Pt_tm1 for the next loop
//
   zt_tm1 = zt_tm1 + (Pt_tm1/(sig2+phi*Pt_tm1*phi'))*phi'*(ut(ti1-1)-phi*zt_tm1);
   Pt_tm1 = Pt_tm1 - Pt_tm1*phi'*phi*Pt_tm1/(sig2+phi*Pt_tm1*phi') + V;
end

disp(' ')
disp('Computing time (seconds) for this task:')
compute_seconds = toc

figure
plot([pcet(maxnlags+1:end)  xtm1])
title('Inflation: actual versus one-step prediction')

figure
plot(pcet(maxnlags+1:end)-xtm1)
title('Inflation: one-step forecast errors')

figure
plot([ut(maxnlags+1:end)  utpred])
title('Unemployment: actual versus one-step prediction')

figure
plot(ut(maxnlags+1:end)-utpred)
title('Unemployment: one-step forecast errors')

jnk1=pcet(maxnlags+1:end)-xtm1;
jnk2=ut(maxnlags+1:end)-utpred;
jnk = [jnk2 jnk1];
save swz3.prn jnk  -ascii