(* ::Package:: *)

(* eztpmodelbase6zc.m

  This Mathematica program file computes the macroeconomic and finance moments for the paper, "The Bond Premium in a
   DSGE Model with Long-Run Real and Nominal Risks," by Glenn Rudebusch and Eric Swanson, forthcoming in the American
   Economic Journal: Macroeconomics.

  First, start a Mathematica session and load Perturbation AIM.  Then run this file.  This file defines the model; calls 
   Perturbation AIM to compute the model's steady state, first-order, second-order, and third-order solutions; uses the
   computed stochastic discount factor solution to solve for zero-coupon bond prices out to 40 quarters; and finally,
   computes the macroeconomic and finance moments of the model by a third-order approximation method around the
   nonstochastic steady state.  Additional details regarding each of these steps is given below.

  This file has been tested in Mathematica 7 on Windows XP and Linux platforms.  It should also work in Mathematica 8 and
   Mathematica 6, although it has not been tested there.

  Written by Eric Swanson, April 2008 to August 2010.
*)

(* First, define the equations of the model, just as in Appendix B of the paper.  Note that the Epstein-Zin objective has the form
     V[t] = C[t]^(1-phi) /(1-phi) + chi0 *(LMax - L[t])^(1-chi) /(1-chi) + beta *(E_t V[t+1]^(1-alpha))^(1/(1-alpha)), which must
     be defined in three equations as described below.  The variable Aperm denotes the growth rate of the technology trend (called Z
     in the paper).  All the variables in the model have already been "normalized" by the (growing) level of technology, so that
     the model as written has a nonstochastic steady state.
*)

eqns = {
(* Epstein-Zin Value Function and Consumer's Euler Equation *)
V[t] == C[t]^(1-phi) /(1-phi) + chi0 *(LMax - L[t])^(1-chi) /(1-chi) + beta *Vkp[t],
C[t]^-phi == beta *(Exp[Int[t]]/pi[t+1]) *C[t+1]^-phi *Aperm[t+1]^-phi *(V[t+1]/Vkp[t] *Aperm[t+1]^(1-phi))^-alphakp, (* Euler eqn *)

(* The following two equations define the E-Z-W-K-P certainty equivalent term Vkp = (E_t V[t+1]^(1-alphakp))^(1/(1-alphakp)).
     It takes two equations to do this because Perturbation AIM sets the expected value of all equations equal to zero,
     E_t F(variables) = 0.  Thus, the first equation below defines Valphaexp[t] == E_t V[t+1]^(1-alphakp).  The second
     equation then takes the (1-alphakp)th root of this expectation.
   Finally, the scaling and unscaling of Valphaexp[t] by the constant VAIMSS improves the numerical behavior of the model; without
     it, the steady-state value of Valphaexp can be minuscule (e.g., 10^-50), which requires Mathematica to use astronomical levels
     of precision to solve.
*)
Valphaexp[t] == (V[t+1]/VAIMSS *Aperm[t+1]^(1-phi)/AgrowBar^(1-phi))^(1-alphakp),
Vkp[t] == VAIMSS*AgrowBar^(1-phi) *Valphaexp[t]^(1/(1-alphakp)),

(* Price-setting equations *)
zn[t] == (1+theta) *MC[t] *Y[t] + xi *beta *(C[t+1]/C[t])^-phi *Aperm[t+1]^-phi *(V[t+1]/Vkp[t] *Aperm[t+1]^(1-phi))^-alphakp
								*pi[t+1]^((1+theta)/theta/eta) *zn[t+1],
zd[t] == Y[t] + xi *beta *(C[t+1]/C[t])^-phi *Aperm[t+1]^-phi *(V[t+1]/Vkp[t] *Aperm[t+1]^(1-phi))^-alphakp *pi[t+1]^(1/theta) *zd[t+1],
p0[t]^(1+(1+theta)/theta *(1-eta)/eta) == zn[t] /zd[t],
pi[t]^(-1/theta) == (1-xi) *(p0[t]*pi[t])^(-1/theta) + xi,

(* Marginal cost and real wage *)
MC[t] == wreal[t] /eta *Y[t]^((1-eta)/eta) /A[t]^(1/eta) /KBar^((1-eta)/eta), (* the extra etas come from industry-specific capital *)
chi0 *(LMax-L[t])^-chi /C[t]^-phi == wreal[t],

(* Output equations *)
Y[t] == A[t] *KBar^(1-eta) *L[t]^eta /Disp[t],
Disp[t]^(1/eta) == (1-xi) *p0[t]^(-(1+theta)/theta/eta) + xi *pi[t]^((1+theta)/theta/eta) *Disp[t-1]^(1/eta),
C[t] == Y[t] - G[t] - IBar,                                               (* aggregate resource constraint, no adj costs for now *)

(* Monetary Policy Rule *)
Log[piavg[t]] == rhoinflavg *Log[piavg[t-1]] + (1-rhoinflavg) *Log[pi[t]],
4*Int[t] == (1-taylrho) * ( 4*Log[1/beta *AgrowBar^phi] + 4*Log[piavg[t]]
	+ taylpi * (4*Log[piavg[t]] - pistar[t]) + tayly * (Y[t]-YBar)/YBar )   (* Taylor rule responds to output gap *)
	+ taylrho * 4*Int[t-1] + epsifl *eps[Int][t],                           (* multiply Int, infl by 4 to put at annual rate *)

(* Exogenous Shocks *)
Log[A[t]] == rhoa * Log[A[t-1]] + epsafl*eps[A][t],
Log[Aperm[t]/AgrowBar] == rhoaperm * Log[Aperm[t-1]/AgrowBar] + epsapermfl*eps[Aperm][t],
Log[G[t]/GBar] == rhog * Log[G[t-1]/GBar] + epsgfl*eps[G][t],
pistar[t] == (1-rhopistar) *piBar + rhopistar *pistar[t-1] + gssload *(4*Log[piavg[t]] - pistar[t]) + epspistarfl*eps[pistar][t],

(* Stochastic Discount Factor (nominal) *)
Intr[t] == Log[Exp[Int[t-1]]/pi[t]], (* ex post real short rate *)
sdf[t] == beta *(C[t]/C[t-1])^-phi *Aperm[t]^-phi *(V[t]/Vkp[t-1] *Aperm[t]^(1-phi))^-alphakp /pi[t]
}


(* These define the parameter values of the model, and are currently set to their "best fit" values.  Thus, this file should
  reproduce the results in the last column of Table 2 of the paper. *)

parametervalues = {
  eta -> 2/3,
  beta -> .99 *AgrowBar^phi,
  delta -> .02,

  IES -> 0.11,
  Frisch -> 0.28,
  CRRA -> 110,

  phi -> 1/IES,
  chi0 -> Exp[wrealAIMSS -phi*CAIMSS] *(LMax-1)^chi, (* normalize L in ss = 1 *)
  chi -> 1/Frisch *(LMax-1),
  alphakp -> (CRRA - 1/(1/phi + 1/chi *(LMax-1))) *(1/(1-phi) + 1/(1-chi) *(LMax-1)),

  theta -> .2,
  xi -> .78,

  AgrowBar -> 1.0025,
  GBar -> .17 *YBar,
  KBar -> 10 *YBar,
  IBar -> (delta +(AgrowBar-1)) *KBar,
  YBar -> Exp[YAIMSS],
  piBar -> 0,
  LMax -> 3,

  taylrho -> .73,
  taylpi -> .53,
  tayly -> .93,

  rhoa -> .96,
  rhoaperm -> 0,
  rhog -> .95,
  rhopistar -> 0,
  epsafl -> 1,    (* these flags turn on or off each corresponding shock in the model *)
  epsapermfl -> 0,
  epsgfl -> 1,
  epsifl -> 1,
  epspistarfl -> 0,
  rhoinflavg -> .7,  (* for price level targeting, set rhoinflavg -> .99 or .999 *)
  gssload -> 0    (* this is the GSS theta parameter in the long-run nominal risks section of the paper *)
}

(* The following lines replace each listed variable X with E^X, so each becomes the log of what it was before *)
loglinearizevars = {A, Aperm, C, Disp, G, L, p0, pi, piavg, sdf, wreal, Y, zd, zn} ;
logRules = Map[#[x_]->E^(#[x])&, loglinearizevars] ;

(* The following two lines define the number of digits of numerical precision to use and the size of a "zero" in the model.
  Given the very high risk aversion and hence curvature of the model, it is wise to use very high numerical precision here, 
  since the numerical errors compound with each additional order of approximation.
*)
AIMPrecision = 90  ; (* 90 digits seems to work well in general, but feel free to use more precision or play around with less. *)
AIMZeroTol = 10^-20 ; (* This says numbers less than 10^-20 should be regarded as being economically zero *)

eztpmodel = eqns /.logRules //.SetPrecision[parametervalues,AIMPrecision] ; (* plug the parameter values into the model *)

(* These are the moments of the shocks, they will be used later *)
momSubs = {Sigma->1, mom[_,3]->0, mom[Aperm,2]->.001^2, mom[A,2]->.005^2, mom[G,2]->.004^2, mom[Int,2]->.003^2, mom[pistar,2]->.0005^2} ;


(* Now, solve for the model's nonstochastic steady state *)

(* A list of all the variables in the model, in alphabetical order, compare to SSGuess:
{A, Aperm, C, Disp, G, Int, Intr, L, MC, p0, pi, piavg, pistar, sdf, V, Valphaexp, Vkp, wreal, Y, zd, zn}
*)
SSGuess = {0, .0025, .65, 0, -.6, .01, .01, 0, .8, 0, 0, 0, 0, 0, -1, 1, -1, .6, 1.15, 2.5, 2.5} ;
ss = AIMSS[eztpmodel,{},AIMSSGuess->SSGuess] ; 
Print[ss];


(* The following lines define some useful abbreviations for variables that we will be using repeatedly below *)
vars = AIMStateVars[eztpmodel] ;
args = AIMGenericArgs[eztpmodel] ;
lagvars = AIMLagVars[eztpmodel] ;
lagvarpos = Flatten[Map[Position[AIMVarNames[eztpmodel],#]&, Map[Head,lagvars]]] ;

sdfpos =  First[Flatten[Position[AIMVarNames[eztpmodel], sdf]]] ;

(* The call to AIMSeries[.,.,3] solves the model to first, second, and then third order *)
sdfsoln3 = AIMSeries[eztpmodel,{},3][[sdfpos,2]] /.x_Real:>N[x] ;

(* Note that sdfsoln3 is for the sdf in logs.  We will need to convert this to levels for bond pricing below. *)
sdfcoeffs = CoefficientArrays[sdfsoln3, vars] ;

sdfargs = Expand[Normal[sdfcoeffs[[1]] + sdfcoeffs[[2]].args + sdfcoeffs[[3]].args.args + sdfcoeffs[[4]].args.args.args]] ;

sdfss = sdfcoeffs[[1]] ; (* steady-state log sdf *)
sdflevel = Exp[sdfss] *Expand[ 1 + (sdfargs-sdfss) + (sdfargs-sdfss)^2 /2 + (sdfargs-sdfss)^3 /6] ;
sdfcoeffs = CoefficientArrays[sdflevel, args] ;

(* Abbreviations for the first-, second-, and third-order solutions to the model, use standard precision now that the model
  has been solved successfully *)
soln1 = N[AIMSoln[eztpmodel,{},1,AIMZeroTol,AIMPrecision]] ;
soln2 = AIMSoln[eztpmodel,{},2,AIMZeroTol,AIMPrecision] /.x_Real:>N[x] ;
soln3 = AIMSoln[eztpmodel,{},3,AIMZeroTol,AIMPrecision] /.x_Real:>N[x] ;

lagsoln1 = soln1[[lagvarpos]] ;
lagsoln2 = soln2[[lagvarpos]] ;
lagsoln3 = soln3[[lagvarpos]] ;


(* Using the sdf solution from the model above, we can now compute the zero-coupon bond price for any maturity bond, by
  iterating on the relationship p_t^(n) = E_t m_{t+1} p_{t+1}^(n-1).  In Mathematica syntax, this is essentially:
    pb2 = Expand[sdfsoln3 *pb1 /.Thread[lagvars->lagsoln3]] /.eps[x_][t]^2->Sigma^2 mom[x,2] /.eps[_][_]->0 /.Sigma->Last[args] ;
  The following function is uglier, but computes the 3rd-order bond price much more quickly.  The function first computes
  the constant term, then the first-order terms (times vars) with the 3rd-order solution substituted in, then the
  second-order terms (times vars.vars) with the 2nd-order solution substituted in, then the third-order terms (times
  vars.vars.vars) with the 1st-order solution substituted in.  This produces a solution that includes all terms up to
  third order, plus some extraneous 4th-order terms, but those last 4th-order terms can easily be dropped by excluding
  the last term of the CoefficientArrays[] result.  Finally, the conditional expectation is taken.
*)

pbfunction[coeffs_] := CoefficientArrays[
        Chop[Expand[sdfcoeffs[[1]]*coeffs[[1]] +
               ((sdfcoeffs[[1]]*coeffs[[2]] + sdfcoeffs[[2]]*coeffs[[1]]).vars /.Thread[lagvars->lagsoln3]) +
               ((sdfcoeffs[[1]]*coeffs[[3]] + Outer[Times,Normal[sdfcoeffs[[2]]],Normal[coeffs[[2]]]]
                   + sdfcoeffs[[3]]*coeffs[[1]]).vars.vars /.Thread[lagvars->lagsoln2]) +
               ((sdfcoeffs[[1]]*coeffs[[4]] + Outer[Times,Normal[sdfcoeffs[[2]]],Normal[coeffs[[3]]]] +
                   Outer[Times,Normal[sdfcoeffs[[3]]],Normal[coeffs[[2]]]] + sdfcoeffs[[4]]*coeffs[[1]]).vars.vars.vars /.
                   Thread[lagvars->lagsoln1])]]
               /.eps[x_][t]^2->Sigma^2 mom[x,2] /.eps[_][_]->0 /.mom[_,3]->0 /.Sigma->Last[args],
        args] ;

nbond = 40 ; (* longest-maturity bond to consider *)
pbacoeffs = Nest[pbfunction,CoefficientArrays[1 + Last[args]^4, args], nbond-1] ; (* nbond-1 period zero-coupon bond price, to compute ehpr *)
pbcoeffs = pbfunction[pbacoeffs] ; (* nbond period zero-coupon bond price *)
pbargs = Expand[Normal[pbcoeffs[[1]] + pbcoeffs[[2]].args + pbcoeffs[[3]].args.args + pbcoeffs[[4]].args.args.args]] ;

(* Note, the initial value in the Nest[] should just be CoefficientArrays[1], but CoefficientArrays[1] is a list with
  only one element (the constant term), instead of the usual four or more elements corresponding to terms from 0th to
  3rd or higher order.  CoefficientArrays[1 + Last[args]^4] solves this problem by adding an extraneous 4th-order
  term that gets dropped anyway.
*)

pbss = sdfcoeffs[[1]]^nbond ; (* steady-state bond price *)
intn = -1/nbond *Expand[ Log[pbss] + 1/pbss *(pbargs-pbss) -1/pbss^2 *(pbargs-pbss)^2 /2 +1/pbss^3 *(pbargs-pbss)^3 /3] ;
intncoeffs = CoefficientArrays[intn, args] ;
intn = Expand[Normal[intncoeffs[[1]] + intncoeffs[[2]].vars + intncoeffs[[3]].vars.vars + intncoeffs[[4]].vars.vars.vars]] ;


(* To compute the term premium, we also need the risk-neutral bond price, computed recursively as p_t^(n) = E_t m_{t+1} E_t p_{t+1}^(n-1). *)

pb1coeffs = pbfunction[CoefficientArrays[1 + Last[args]^4, args]] ;

pbexpfunction[coeffs_] := CoefficientArrays[
        Chop[Expand[coeffs[[1]] +
               (coeffs[[2]].vars /.Thread[lagvars->lagsoln3]) +
               (coeffs[[3]].vars.vars /.Thread[lagvars->lagsoln2]) +
               (coeffs[[4]].vars.vars.vars /.Thread[lagvars->lagsoln1])]]
               /.eps[x_][t]^2->Sigma^2 mom[x,2] /.eps[_][_]->0 /.mom[_,3]->0 /.Sigma->Last[args],
        args] ;

pbrnfunction[coeffs_] := CoefficientArrays[pbexpcoeffs = pbexpfunction[coeffs] ;
        Chop[Expand[pb1coeffs[[1]]*pbexpcoeffs[[1]] +
               (pb1coeffs[[1]]*pbexpcoeffs[[2]] + pb1coeffs[[2]]*pbexpcoeffs[[1]]) .args +
               (pb1coeffs[[1]]*pbexpcoeffs[[3]] + Outer[Times,Normal[pb1coeffs[[2]]],Normal[pbexpcoeffs[[2]]]]
                   + pb1coeffs[[3]]*pbexpcoeffs[[1]]).args.args +
               (pb1coeffs[[1]]*pbexpcoeffs[[4]] + Outer[Times,Normal[pb1coeffs[[2]]],Normal[pbexpcoeffs[[3]]]]
                   + Outer[Times,Normal[pb1coeffs[[3]]],Normal[pbexpcoeffs[[2]]]] + pb1coeffs[[4]]*pbexpcoeffs[[1]]).args.args.args]],
        args] ; 

pbrncoeffs = Nest[pbrnfunction,pb1coeffs,nbond-1] ;
pbrnargs = Expand[Normal[pbrncoeffs[[1]] + pbrncoeffs[[2]].args + pbrncoeffs[[3]].args.args + pbrncoeffs[[4]].args.args.args]] ;

intrn = -1/nbond *Expand[ Log[pbss] + 1/pbss *(pbrnargs-pbss) -1/pbss^2 *(pbrnargs-pbss)^2 /2 +1/pbss^3 *(pbrnargs-pbss)^3 /3] ;
intrncoeffs = CoefficientArrays[intrn, args] ;
intrn = Expand[Normal[intrncoeffs[[1]] + intrncoeffs[[2]].vars + intrncoeffs[[3]].vars.vars + intrncoeffs[[4]].vars.vars.vars]] ;


(* print out the steady-state term premium estimate: *)
tp = 40000 *Chop[intn - intrn] ;
meantp = tp /.momSubs /.Thread[lagvars->0] /.eps[_][_]->0 ;
Print["meantp = ", meantp ] ;


(* Now compute the macro and finance moments of the model to third order.  One can do this in principle by simulation, but the
     simulations sometimes explode, especially for more curved models.  The following blocks of code compute the macro and
     finance moments of the model to third order, without any need for simulation. *)

(* Compute the AIM Bmat and S0Inv matrices.  Bmat is the first-order dynamics of the variables, and S0Inv is the
     contemporaneous response of the variables to shocks *)
eqnseq0 = Flatten[{eztpmodel /.x_Real:>SetPrecision[x,AIMPrecision]}] /.Equal->Subtract ;
allvars = Flatten[Map[Through[AIMVarNames[eztpmodel] [t+#]]&, Range[-Max[AIMMaxLag[eztpmodel],1], AIMMaxLead[eztpmodel]]]] ;
hmat = Chop[Outer[D,eqnseq0,allvars] /.AIMSSSubs /.ss, AIMZeroTol] ;
epsmat = Chop[-Outer[D,eqnseq0,AIMShocks[eztpmodel]] /.AIMSSSubs /.ss,AIMZeroTol] ;
{cofb,s0inv} = AIMLinearSoln[hmat,AIMMaxLead[eztpmodel]] ;
s0invsmall = N[s0inv .epsmat] ;


(* Start computing unconditional moments *)
momchop = 10^-7 ;

(* First, define some useful placekeeping variables *)
xpos = Range[Length[lagvars]] ;
epspos = Range[Length[lagvars]+1,Length[AIMStateVars[eztpmodel]]] ;
varsx = ReplacePart[AIMStateVars[eztpmodel],Thread[epspos->0]] ;
varseps = ReplacePart[AIMStateVars[eztpmodel],Thread[xpos->0]] ;


(* Compute Ex^2 to 2nd order *)
Print["Computing Ex^2 to 2nd order, time is ", Date[]] ;
sigmaeps = s0invsmall .DiagonalMatrix[AIMShocks[eztpmodel] /.eps[x_][_]->mom[x,2] /.momSubs] .Transpose[s0invsmall] ;
b = N[cofb] ;
sigmax = FixedPoint[(s = # + b.#.Transpose[b]; b = b.b; s)&, sigmaeps, 8] ; (* doubling algorithm *)
subsorder2 = Intersection[Thread[Flatten[Outer[Times,lagvars,lagvars]] -> Flatten[sigmax[[lagvarpos,lagvarpos]]]]] ;


(* Compute Ex to 2nd order *)
Print["Computing Ex to 2nd order, time is ", Date[]] ;
mua = soln2 /.momSubs ;
mucoeffs = CoefficientArrays[mua,args] ;

mub = Expand[Normal[mucoeffs[[3]] .vars.vars]] /.eps[x_][t]^2->mom[x,2] /.eps[_][_]->0 /.momSubs /.subsorder2 ;

b = N[cofb] ;
mu = FixedPoint[(s = # + b.# ; b = b.b; s)&, mub, 8] ; (* doubling *)  (* NOTE: need to add AIMSS if you want levels! *)

subsorder1 = Thread[lagvars->mu[[lagvarpos]]] ;


(* Compute E x^4 to 4th order *)
Print["Computing Ex^4 to 4th order, time is ", Date[]] ;
xfoura = Expand[Outer[Times,lagsoln1,lagsoln1,lagsoln1,lagsoln1]] ;
xfourcoeffs = CoefficientArrays[xfoura,args[[xpos]]] ;

xfourb = Expand[Normal[xfourcoeffs[[3]].lagvars.lagvars]] /.Thread[args[[epspos]]->vars[[epspos]]] /.
					eps[x_][t]^2->mom[x,2] /.eps[_][_]->0 /.momSubs /.subsorder2 ;
xfourc = Expand[Normal[xfourcoeffs[[1]]]] /.Thread[args[[epspos]]->vars[[epspos]]] /.
					eps[x_][t]^4->3*mom[x,2]^2 /.eps[x_][t]^2->mom[x,2] /.eps[_][_]->0 /.momSubs ;

b = N[cofb][[lagvarpos,lagvarpos]] ;
xfour = FixedPoint[(s = # + Transpose[b .Transpose[b .Transpose[b. Transpose[b. #,
							{2,3,4,1}], {2,3,4,1}], {2,3,4,1}], {2,3,4,1}]; b = b.b; s)&,
			xfourb + xfourc, 8] ; (* doubling algorithm for 4-dimensional tensor *)
subsorder4 = Intersection[Thread[Flatten[Outer[Times,lagvars,lagvars,lagvars,lagvars]] -> Flatten[xfour]] ] ;


(* Compute E x^3 to 4th order *)
Print["Computing Ex^3 to 4th order, time is ", Date[]] ;
xsoln2 = lagsoln2 /.momSubs ;
xthreea = Expand[Outer[Times,xsoln2,xsoln2,xsoln2]] ;
xthreecoeffs = CoefficientArrays[xthreea,args] ;

xthreeb = Chop[Expand[Normal[xthreecoeffs[[4]].vars.vars.vars]] - Expand[Normal[xthreecoeffs[[4]].varsx.varsx.varsx]]] /.
					eps[x_][t]^2->mom[x,2] /.eps[_][_]->0 /.momSubs /.subsorder1 ;

xthreec = Expand[Normal[xthreecoeffs[[5]].vars.vars.vars.vars]] /.
					eps[x_][t]^4->3*mom[x,2]^2 /.eps[x_][t]^2->mom[x,2] /.eps[_][_]->0 /.momSubs /.
					subsorder4 /.subsorder2 ;

b = N[cofb[[lagvarpos,lagvarpos]]] ;
xthree = FixedPoint[(s = # + Transpose[b .Transpose[b .Transpose[b. #, {2,3,1}], {2,3,1}], {2,3,1}]; b = b.b; s)&,
			xthreeb + xthreec, 8] ; (* doubling algorithm for 3-dimensional tensor *)
subsorder3 = Intersection[Thread[Flatten[Outer[Times,lagvars,lagvars,lagvars]] -> Flatten[xthree]] ] ;


(* Solve for E x^2 to 4th order *)
reportvars = Union[Map[Head,AIMLagVars[eztpmodel]], {C, Y, L, wreal, pi, Int, Intr}] ;
reportvarpos = Flatten[Map[Position[AIMVarNames[eztpmodel],#]&, reportvars]] ;
lagvarposa = Flatten[Map[Position[reportvars,#]&, Map[Head,lagvars]]] ;
Print["Solving for Ex^2 to 4th order, time is ", Date[]] ;
xsoln3 = Chop[soln3[[reportvarpos]], momchop] /.momSubs ;
xtwoa = Expand[Outer[Times,xsoln3,xsoln3]] ;
xtwocoeffs = CoefficientArrays[xtwoa,args] ;

xtwob =  Expand[Normal[xtwocoeffs[[3]].varseps.varseps]] /.eps[x_][t]^2->mom[x,2] /.eps[_][_]->0 /.momSubs ;

xtwoc =  Expand[Normal[xtwocoeffs[[4]].vars.vars.vars]] /.eps[x_][t]^2->mom[x,2] /.eps[_][_]->0 /.momSubs /.subsorder3 /.subsorder1 ;

xtwod =  Expand[Normal[xtwocoeffs[[5]].vars.vars.vars.vars]] /.eps[x_][t]^4->3*mom[x,2]^2 /.eps[x_][t]^2->mom[x,2] /.eps[_][_]->0 /.
					momSubs /.subsorder4 /.subsorder2 ;

b = N[cofb][[reportvarpos,reportvarpos]] ;
sigmaxfour = FixedPoint[(s = # + b.#.Transpose[b]; b = b.b; s)&, xtwob + xtwoc + xtwod, 8] ; (* doubling *)
subsorder2a = Intersection[Thread[Flatten[Outer[Times,lagvars,lagvars]] -> Flatten[sigmaxfour[[lagvarposa,lagvarposa]]]]] ;
sigmaxfour = sigmaxfour - Outer[Times,mu[[reportvarpos]],mu[[reportvarpos]]] ; (* subtract means *)


(* Print moments to 2nd and 4th order :*)
Print["Theoretical std deviations of {C, Y, L, wreal, pi, Int, Intr} to 1st order:"] ;
moments = Sqrt[Tr[sigmax,List]] [[
	Flatten[Map[Position[AIMVarNames[eztpmodel], #] &, {C, Y, L, wreal, pi, Int, Intr}]] ]
	] * {100, 100, 100, 100, 400, 400, 400}
Print[moments] ;

Print["Theoretical std deviations of {C, Y, L, wreal, pi, Int, Intr} to 2nd order:"] ;
moments4 = Sqrt[Tr[sigmaxfour,List]] [[
	Flatten[Map[Position[reportvars, #] &, {C, Y, L, wreal, pi, Int, Intr}]] ]
	] * {100, 100, 100, 100, 400, 400, 400}
Print[moments4] ;


(* Compute mean of term premium to 3rd order and variance of term premium to 6th order *)
Etp2 = Expand[(tp /.momSubs)^2] /.eps[x_][t]^2->mom[x,2] /.eps[_][_]->0 /.momSubs /.subsorder2 /.subsorder1 ;
stdtp = Sqrt[Etp2 - meantp^2] ; (* subtract mean *)


(* Compute mean of long-term bond yield to 2nd order and variance of long-term bond yield to 4th order *)
intnord2 = Expand[intncoeffs[[1]] + intncoeffs[[2]].vars + intncoeffs[[3]].vars.vars] ;        (* intn to 2nd order *)
Eintn = intnord2 /.eps[x_][t]^2->mom[x,2] /.eps[_][_]->0 /.momSubs /.subsorder2 /.subsorder1 ; (* Eintn to 2nd order *)
intnss = intncoeffs[[1]] ;
intn2coeffs = CoefficientArrays[(intn-intnss)^2,vars] ;
Eintn2 = (Expand[intn2coeffs[[3]].vars.vars] /.eps[x_][t]^2->mom[x,2] /.eps[_][_]->0 /.momSubs /.subsorder2a) +
            (Expand[intn2coeffs[[4]].vars.vars.vars] /.eps[x_][t]^2->mom[x,2] /.eps[_][_]->0 /.momSubs /.subsorder3 /.subsorder1) + 
            (Expand[intn2coeffs[[5]].vars.vars.vars.vars] /.                   (* Terms in Eintn2 are all accurate to 4th order *)
                                eps[x_][t]^4->3*mom[x,2]^2 /.eps[x_][t]^2->mom[x,2] /.eps[_][_]->0 /.momSubs /.subsorder4 /.subsorder2) ;
stdintn = Sqrt[Eintn2 - (Eintn-intnss)^2] ; (* Eintn-intnss =0 to 1st order, so (Eintn-intnss)^2 is accurate to 4th order *)


(* Compute mean of slope to 2nd order and variance of slope to 4th order *)
intpos = First[Flatten[Position[AIMVarNames[eztpmodel],Int]]] ;
slopeord2 = Chop[intnord2 - (AIMSeries[eztpmodel,2][[intpos,2]] /.x_Real:>N[x])] ;               (* slope to 2nd order *)
Eslope = slopeord2 /.eps[x_][t]^2->mom[x,2] /.eps[_][_]->0 /.momSubs /.subsorder2 /.subsorder1 ; (* Eslope to 2nd order *)
slope = Chop[intn - (AIMSeries[eztpmodel,3][[intpos,2]] /.x_Real:>N[x])] ;                       (* slope to 3rd order *)
slope2coeffs = CoefficientArrays[slope^2,vars] ;
Eslope2 = (Expand[slope2coeffs[[3]].vars.vars] /.eps[x_][t]^2->mom[x,2] /.eps[_][_]->0 /.momSubs /.subsorder2a) +
            (Expand[slope2coeffs[[4]].vars.vars.vars] /.eps[x_][t]^2->mom[x,2] /.eps[_][_]->0 /.momSubs /.subsorder3 /.subsorder1) + 
            (Expand[slope2coeffs[[5]].vars.vars.vars.vars] /.                  (* Terms in Eslope2 are all accurate to 4th order *)
                                eps[x_][t]^4->3*mom[x,2]^2 /.eps[x_][t]^2->mom[x,2] /.eps[_][_]->0 /.momSubs /.subsorder4 /.subsorder2) ;
stdslope = Sqrt[Eslope2 - Eslope^2] ; (* Eslope = 0 to 1st order, so Eslope^2 is accurate to 4th order *)


(* Compute mean excess holding period return to 2nd order and variance of excess holding period return to 4th order.
  Thus, we want to compute 3rd-order expression for:  ehpr_t = pb_t^(39)/pb_{t-1}^(40) - Exp[Int[t-1]]
                                                             = pb_t^(39)/pb_{t-1}^(40) - 1/pb_{t-1}^(1) ;
  Since we are interested in unconditional moments, we can work with E ehpr_{t+1} and Var ehpr_{t+1} instead.  This produces
  expressions with eps[_][t+1] as well as eps[_][t] terms, which we must keep track of in computing the mean and variance.
*)
pbinv = 1/pbss -1/pbss^2 *(pbargs-pbss) + 1/pbss^3 *(pbargs-pbss)^2 -1/pbss^4 *(pbargs-pbss)^3 ;
pbinvcoeffs = CoefficientArrays[pbinv,args] ;
pbinvord2 = Expand[pbinvcoeffs[[1]] + pbinvcoeffs[[2]].vars + pbinvcoeffs[[3]].vars.vars] ;
pbinv = pbinvord2 + Expand[pbinvcoeffs[[4]].vars.vars.vars] ;

pb1ss = pb1coeffs[[1]] ;
pb1args = Expand[pb1coeffs[[1]] + pb1coeffs[[2]].args + pb1coeffs[[3]].args.args + pb1coeffs[[4]].args.args.args] ;
pb1inv = 1/pb1ss -1/pb1ss^2 *(pb1args-pb1ss) + 1/pb1ss^3 *(pb1args-pb1ss)^2 -1/pb1ss^4 *(pb1args-pb1ss)^3 ;
pb1invcoeffs = CoefficientArrays[pb1inv,args] ;
pb1invord2 = Expand[pb1invcoeffs[[1]] + pb1invcoeffs[[2]].vars + pb1invcoeffs[[3]].vars.vars] ;
pb1inv = pb1invord2 + Expand[pb1invcoeffs[[4]].vars.vars.vars] ;

varstpo1 = vars /.Thread[lagvars->lagsoln1] /.t->t+1 /.Thread[args->vars] ;
varstpo2 = vars /.Thread[lagvars->lagsoln2] /.t->t+1 /.Thread[args->vars] ;
varstpo3 = vars /.Thread[lagvars->lagsoln3] /.t->t+1 /.Thread[args->vars] ;

ehprord2 = Chop[Expand[pbacoeffs[[1]]*pbinvcoeffs[[1]] +
               pbacoeffs[[1]]*pbinvcoeffs[[2]].vars + pbacoeffs[[2]].varstpo2*pbinvcoeffs[[1]] +
               pbacoeffs[[1]]*pbinvcoeffs[[3]].vars.vars + pbacoeffs[[2]].varstpo1 *pbinvcoeffs[[2]].vars
                   + pbacoeffs[[3]].varstpo1.varstpo1*pbinvcoeffs[[1]]] - pb1invord2] ;

Eehpr = ehprord2 /.eps[x_][_]^2->Sigma^2 mom[x,2] /.eps[_][_]->0 /.momSubs /.subsorder2 /.subsorder1 ; (* E ehpr to 2nd order *)

ehpr = Chop[Expand[pbacoeffs[[1]]*pbinvcoeffs[[1]] +
               pbacoeffs[[1]]*pbinvcoeffs[[2]].vars + pbacoeffs[[2]].varstpo3*pbinvcoeffs[[1]] +
               pbacoeffs[[1]]*pbinvcoeffs[[3]].vars.vars + pbacoeffs[[2]].varstpo2 *pbinvcoeffs[[2]].vars
                   + pbacoeffs[[3]].varstpo2.varstpo2*pbinvcoeffs[[1]] +
               pbacoeffs[[1]]*pbinvcoeffs[[4]].vars.vars.vars
                   + pbacoeffs[[2]].varstpo1 *pbinvcoeffs[[3]].vars.vars
                   + pbacoeffs[[3]].varstpo1.varstpo1 *pbinvcoeffs[[2]].vars
                   + pbacoeffs[[4]].varstpo1.varstpo1.varstpo1 *pbinvcoeffs[[1]]] - pb1inv] ;

ehprvars = Join[vars,AIMShocks[eztpmodel]/.t->t+1] ;
ehpr2coeffs = CoefficientArrays[ehpr^2,ehprvars] ;
Eehpr2 = (Expand[ehpr2coeffs[[3]].ehprvars.ehprvars] /.eps[x_][_]^2->mom[x,2] /.eps[_][_]->0 /.momSubs /.subsorder2a) +
            (Expand[ehpr2coeffs[[4]].ehprvars.ehprvars.ehprvars] /.eps[x_][t]^2->mom[x,2] /.eps[_][_]->0 /.momSubs /.subsorder3 /.subsorder1) +
            (Expand[ehpr2coeffs[[5]].ehprvars.ehprvars.ehprvars.ehprvars] /.        (* Terms in Eehpr2 are all accurate to 4th order *)
                                eps[x_][_]^4->3*mom[x,2]^2 /.eps[x_][_]^2->mom[x,2] /.eps[_][_]->0 /.momSubs /.subsorder4 /.subsorder2) ;
stdehpr = Sqrt[Eehpr2 - Eehpr^2] ; (* Eehpr = 0 to 1st order, so Eehpr^2 is accurate to 4th order *)



Print["Theoretical std deviations of {C, Y, L, wreal, pi, Int, Intr, intn} to 2nd order, w/ sd of tp to 3rd order:"] ;
moments46 = Join[Sqrt[Tr[sigmaxfour,List]] [[
	Flatten[Map[Position[reportvars, #] &, {C, Y, L, wreal, pi, Int, Intr}]] ]
	] * {100, 100, 100, 100, 400, 400, 400}, {stdintn}*400, {stdtp}] ;
Print[moments46] ;

Print["mean {termprem, slope, ehpr} to 2nd order (in bp): ", meantp,", ", Eslope*40000,", ", Eehpr*40000] ;
Print["std deviation of {termprem, slope, ehpr} to 2nd order (in pct): ", stdtp/100,", ", stdslope*400,", ", stdehpr*400] ;

momentlist = Join[Drop[moments46,-1],{meantp/100,Eslope*400,Eehpr*400,stdtp/100,stdslope*400,stdehpr*400}] ;

momenttargets = {.83, 1.50, 1.71, .82, 2.52, 2.71, 2.30, 2.41, 1.06, 1.43, 1.76, .54, 1.33, 23.43} ;
momentweights = {1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .1^2} ;
macromomentweights = {1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0} ;
financemomentweights = {0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, .1^2} ;
distance = ((momentlist-momenttargets) *momentweights) .(momentlist-momenttargets) ;
macrodistance = ((momentlist-momenttargets) *macromomentweights) .(momentlist-momenttargets) ;
financedistance = ((momentlist-momenttargets) *financemomentweights) .(momentlist-momenttargets) ;



(* Uncomment the following lines to export matrices of first-order macro dynamics and third-order term premium and long-term bond
  price dynamics.  These matrices can then be used to compute impulse responses, as in the Matlab program, figure2.m *)
(* 
Export["figures/eztpmodelcofb.txt", cofb, "Table"] ; (* first-order dynamics of macro variables *)
Export["figures/eztpmodels0inv.txt", s0inv, "Table"] ; (* first-order contemporaneous response of macro variables to shocks *)

allvarslag = Through[AIMVarNames[eztpmodel] [t-1]] ;
tpcofb = Map[D[tp /.momSubs,#]&, allvarslag] /.Thread[vars->0] ;

tps0invsmall = Map[D[tp /.momSubs,#]&, AIMShocks[eztpmodel]] /.Thread[vars->0] ;
tps0inv = tps0invsmall .Transpose[epsmat] ;

Export["figures/eztpmodeltpcofb.txt", tpcofb, "Table"] ; (* third-order dynamics of term premium *)
Export["figures/eztpmodeltps0inv.txt", tps0inv, "Table"] ; (* third-order contermporaneous response of term premium to shocks *)


pbcofb = Map[D[pbargs /.Thread[args->vars] /.momSubs,#]&, allvarslag] /.Thread[lagvars->0] /.Thread[AIMShocks[eztpmodel]->0] ;
pbs0invsmall = Map[D[pbargs /.Thread[args->vars] /.momSubs,#]&, AIMShocks[eztpmodel]] /.Thread[lagvars->0] /.Thread[AIMShocks[eztpmodel]->0] ;
pbs0inv = pbs0invsmall .Transpose[epsmat] ;

Export["figures/eztpmodelpbcofb.txt", pbcofb, "Table"] ; (* first-order dynamics of long-term bond price *)
Export["figures/eztpmodelpbs0inv.txt", pbs0inv, "Table"] ; (* first-order contemporaneous response of long-term bond price to shocks *)
*)
