1// This file replicates the estimation of the CIA model from 2// Frank Schorfheide (2000) "Loss function-based evaluation of DSGE models" 3// Journal of Applied Econometrics, 15, 645-670. 4// the data are the ones provided on Schorfheide's web site with the programs. 5// http://www.econ.upenn.edu/~schorf/programs/dsgesel.ZIP 6// You need to have fsdat.m in the same directory as this file. 7// This file replicates: 8// -the posterior mode as computed by Frank's Gauss programs 9// -the parameter mean posterior estimates reported in the paper 10// -the model probability (harmonic mean) reported in the paper 11// This file was tested with dyn_mat_test_0218.zip 12// the smooth shocks are probably stil buggy 13// 14// The equations are taken from J. Nason and T. Cogley (1994) 15// "Testing the implications of long-run neutrality for monetary business 16// cycle models" Journal of Applied Econometrics, 9, S37-S70. 17// Note that there is an initial minus sign missing in equation (A1), p. S63. 18// 19// Michel Juillard, February 2004 20options_.usePartInfo=0; 21var m P c e W R k d n l gy_obs gp_obs Y_obs P_obs y dA P2 c2; 22varexo e_a e_m; 23 24parameters alp bet gam mst rho psi del; 25 26alp = 0.33; 27bet = 0.99; 28gam = 0.003; 29mst = 1.011; 30rho = 0.7; 31psi = 0.787; 32del = 0.02; 33 34model ; 35dA = exp(gam+e_a); 36log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m; 37-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c2(+1)*P2(+1)*m(+1))=0; 38W = l/n; 39-(psi/(1-psi))*(c*P/(1-n))+l/n = 0; 40R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W; 411/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0; 42c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1); 43P*c = m; 44m-1+d = l; 45e = exp(e_a); 46y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a)); 47gy_obs = dA*y/y(-1); 48gp_obs = (P/P(-1))*m(-1)/dA; 49Y_obs/Y_obs(-1) = gy_obs; 50P_obs/P_obs(-1) = gp_obs; 51P2 = P(+1); 52c2 = c(+1); 53end; 54 55steady_state_model; 56 dA = exp(gam); 57 gst = 1/dA; 58 m = mst; 59 khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1)); 60 xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1); 61 nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp ); 62 n = xist/(nust+xist); 63 P = xist + nust; 64 k = khst*n; 65 66 l = psi*mst*n/( (1-psi)*(1-n) ); 67 c = mst/P; 68 d = l - mst + 1; 69 y = k^alp*n^(1-alp)*gst^alp; 70 R = mst/bet; 71 W = l/n; 72 ist = y-c; 73 q = 1 - d; 74 75 e = 1; 76 77 gp_obs = m/dA; 78 gy_obs = dA; 79 Y_obs = 1; 80 P_obs = 1; 81 P2 = P; 82 c2 = c; 83end; 84 85 86shocks; 87var e_a; stderr 0.014; 88var e_m; stderr 0.005; 89end; 90 91steady(nocheck); 92 93stoch_simul(aim_solver, order=1, irf=0); 94 95benchmark = load('fs2000_b1L1L_results'); 96threshold = 1e-8; 97 98if max(max(abs(benchmark.oo_.dr.ghx-oo_.dr.ghx) > threshold)); 99 error('error in ghx'); 100elseif max(max(abs(benchmark.oo_.dr.ghu-oo_.dr.ghu) > threshold)); 101 error('error in ghy'); 102end; 103