rss2007model.m

(* Builds on the model from Rudebusch and Swanson (2008, FRB St. Louis Economic Review) and FRBSL conference, with the
   following modifications:
     1. put interest rate rule in terms of output gap rather than output growth
     2. modify the consol to give it a duration of 10 yrs instead of 25 yrs
     3. adjust parameter values to get steady-state ratios roughly calibrated to US data
     4. introduce an adjustment cost parameter psi that we'll set to 0 here but play
	around with later.

Eric Swanson, 8/2007.
*)

(*Barton Baker copy *)




(* objective = (C[t] - habitsize *(C[t-1]))^(1-phi) /(1-phi) - chi0 *L[t]^(1+chi) /(1+chi)
 *)

eqns = {

(* Marginal Utility of Consumption and Consumer's Euler Equation *)
MUc[t] == (C[t] - b*C[t-1])^(-gamma),


(*MUc[t] == beta *(Exp[Int[t]]/(pi[t+1]+1)) *MUc[t+1],*)
(*Try alternative:*)
MUc[t] == beta *(Exp[Int[t]]/(pi[t+1])) *MUc[t+1],



(* Price-setting equations *)
zn[t] == (1+theta)*chi0/(1-alpha)*L[t]^(1+chi)*Disp[t]^(-1/(1-alpha))+beta*xi*pi[t+1]^(((1+theta)/theta)*(1/(1-alpha)))*zn[t+1],
(*zn[t] == (1+theta) *MUc[t] *MC[t] *Y[t] + beta *xi *pi[t+1]^((1+theta)/theta/(1-alpha)) *zn[t+1],*)

zd[t] == Y[t]*MUc[t]+beta*xi*pi[t+1]^(1/theta)*zd[t+1],
(*zd[t] == Y[t] *MUc[t] + beta *xi *pi[t+1]^(1/theta) *zd[t+1],*)

(*This is unchanged *)
p0[t]^(1+((1+theta)/theta)*(alpha/(1-alpha))) == zn[t]/zd[t],

(*Check this out!!!!!!!!!!!!!!!*)
(*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*)
(*
pi[t]^(-1/theta) == (1-xi) *(p0[t]*pi[t])^(-1/theta) + xi,
*)
(*Try out this simpler version*)
pi[t]^(-1/theta) == xi + (1-xi)*(p0[t-1]/p0[t])^(-1/theta),



(* OLD: Marginal cost and quadratic adjustment costs to labor, psi/2 *Log[L[t]/L[t-1]]^2 *)
(* Marginal cost (maybe don't need this???)*)
(*Just use relationship between real wage*)

(*
MC[t] == wreal[t] /(1-alpha) *Y[t]^(alpha/(1-alpha)) /A[t]^(1/(1-alpha)) /KBar^(alpha/(1-alpha)),
chi0 *L[t]^chi /MUc[t] == wreal[t] - psi /L[t] *Log[L[t]/L[t-1]]
						+ beta *psi /L[t] *Log[L[t+1]/L[t]] *MUc[t+1]/MUc[t],
*)
(* Output equations *)
Y[t] == Disp[t]^(-1) *A[t] *KBar^alpha *L[t]^(1-alpha),


(*Try this out or check math again later*)
(*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*)

Disp[t]^(1/(1-alpha)) == (1-xi) *p0[t]^(-(1+theta)/theta/(1-alpha))
	+ xi *pi[t]^((1+theta)/theta/(1-alpha)) *Disp[t-1]^(1/(1-alpha)),


(*C[t] == Y[t] - G[t] - IBar - psi/2 *Log[L[t]/L[t-1]]^2, (* aggregate resource constraint w/ adj costs *)*)
C[t] == Y[t] - G[t], (*aggregate resource constraint*)


(* Monetary Policy Rule *)

Int[t] == taylrho*Int[t-1]+(1-taylrho)*(Intss+gy(Y[t]-Y[t-1])+gpi*pi[t])+eps[Int][t],

(*Intr[t] == Log[Exp[Int[t-1]]/pi[t]],*)

(*
piavg[t] == rhoinflavg *piavg[t-1] + (1-rhoinflavg) *pi[t],
*)
(*piavg[t] == (pi[t] + pi[t-1] + pi[t-2] + pi[t-3]) /4,*)
(*
4*Int[t] == (1-taylrho) * ( 4*Log[1/beta] + 4*Log[piavg[t]]
	+ taylpi * (4*Log[piavg[t]] - piBar)
	+ tayly*(Y[t]-YBar)/YBar )
	+ taylrho * 4*Int[t-1] + eps[Int][t],  (* multiply Int, infl by 4 to put at annual rate *)
*)

(*Intr[t] == Log[Exp[Int[t-1]]/pi[t]], (* ex post real short rate *)*)

(* Exogenous Shocks *)

(*Log[A[t]/ABar] == rhoa * Log[A[t-1]/ABar] + eps[A][t],*)
Log[A[t]/ABar] == rhoa * Log[A[t-1]/ABar] + eps[A][t],

(*Log[G[t]/GBar] == rhog * Log[G[t-1]/GBar] + eps[G][t],*)

Log[G[t]/GBar] == rhog * Log[G[t-1]/GBar] + eps[G][t],


(*pistar[t] == (1-rhopistar) *piBar + rhopistar *pistar[t-1] + 0*eps[pistar][t],*)

(*Stochastic Discount factor*)

(* Long-Term Bond Price, Yield, and Term Premium *)



pricebond[t] == 1/100 + beta *consoldelta *MUc[t+1] /MUc[t] /pi[t+1] *pricebond[t+1],
(*pricebond[t] == 1 + beta*MUc[t+1] /MUc[t] /pi[t+1] *pricebond[t+1],*)

pricebondrn[t] == 1/100 + consoldelta *pricebondrn[t+1] *Exp[-Int[t]], (* risk-neutral bond price *)
(*pricebondrn[t] == 1 + *pricebondrn[t+1] *Exp[-Int[t]],*)

ytm[t] == Log[consoldelta*pricebond[t]/(pricebond[t]-1/100)] *400,  (* yield in annualized percent *)
(*ytm[t] == Log[pricebond[t]/(pricebond[t]-1)],*)

ytmrn[t] == Log[consoldelta*pricebondrn[t]/(pricebondrn[t]-1/100)] *400,  (* yield in annualized percent *)
(*ytmrn[t] == Log[pricebondrn[t]/(pricebondrn[t]-1)],*)

termprem[t] == 100 * (ytm[t] - ytmrn[t]) (* term premium in annualized basis points *)

(*termprem[t] == 4 * (ytm[t] - ytmrn[t])*)

}


parametervalues = {
  alpha -> .3,
  beta -> .99,
  theta -> .2,
  xi -> .75,
  gamma -> 2,
  b -> .66,
  chi0 -> (1-b)^(-gamma),
  chi -> 1.5,
  gpi -> 2,
  gy -> .5,
  rhoa -> .9,
  rhog -> .9,
  taylrho -> .7,
  KBar -> 1,
  piss -> 0,
  Intss -> Log[(1+piss)/beta],
  ABar -> 1,
  GBar -> .17*YBar,
  YBar -> Exp[YAIMSS],
  consoldelta -> Exp[ytmAIMSS/400] *(1-1/40)  (* consol depreciation, set to make duration of consol =10 yrs *) 

}


  (*loglinearizevars = {A, C, Disp, G, L, MUc, p0, pi, piavg, pricebond, pricebondrn, wreal, Y, zd, zn} ;*)
loglinearizevars = {A, C, Disp, G, L, MUc, p0, pi, piavg, pricebond, pricebondrn, Y, zd, zn} ;
logRules = Map[#[x_]->E^(#[x])&, loglinearizevars] ;

(*OLD: List of all variables:
{A, C, Disp, G, Int, Intr, L, MC, MUc, p0,
  pi, piavg, pistar, pricebond, pricebondrn, termprem, wreal, Y, ytm, ytmrn, zd, zn}
*)

(*List of all variables:
{A, C, Disp, G, Int, Intr,  L, MUc, p0, pi, pricebond, pricebondrn, termprem, Y, ytm, ytmrn, zd, zn}
*)


SSGuess = {0, .5, 0, -.8, .01, 0, 1.1, 0,
  0, -1, -1, 0, .5, 4, 4, 3.5, 3.5}

AIMPrecision = 250 ; (* 250 *)
AIMZeroTol = 10^-10 (* 10^-50 *) ;  (* for psi = 1000, AIMZeroTol = 10^-8 and Precision 80 seems to work *)

rssmodel = eqns /.logRules //.SetPrecision[parametervalues,AIMPrecision] ;

ss = AIMSS[rssmodel,{},AIMSSGuess->SSGuess] ;

momSubs = {Sigma->1, mom[_,3]->0, mom[A,2]->.01^2, mom[G,2]->.004^2, mom[Int,2]->.004^2}


(* print out the steady-state term premium estimate: *)

tppos = Position[AIMVarNames[rssmodel], termprem][[1,1]] ;

Print[AIMSeries[rssmodel,2][[tppos]]/ /.momSubs /.x_Real:> N[x]] ;



(* Simulate solution forward nperiods periods *)

Print["\nThe Part[] function will issue some warnings here; fear not, everything is ok..."]

nperiods = 1000 ; (* number of periods *)
simdegree = 3 ;  (* degree of approximation *)

reportvars = {C, Y, L, pi, Int, ytm, termprem} ; (* variables of interest *)
(*reportvars = {C, Y, L, wreal, pi, Int, Intr, ytm, termprem} ; (* variables of interest *) *)
simvars = Join[ Map[Head, AIMLagVars[rssmodel]], reportvars, {pricebond, ytmrn}] ; (* variables to simulate *)

lagvarpos = Range[Length[AIMLagVars[rssmodel]]] ;
reportvarpos = Length[AIMLagVars[rssmodel]] + Range[Length[reportvars]] ;

(*shocks = RandomReal[NormalDistribution[0,1],{nperiods,Length[AIMShocks[rssmodel]]}] .
					DiagonalMatrix[AIMShocks[rssmodel] /.eps[x_][_]->Sqrt[mom[x,2]] /.momSubs] ;*)

Here is where we tell it what shock to use

shocks= ConstantArray[0,{nperiods,Length[AIMShocks[rssmodel]]}];

shocks[[1,1]]=.01

soln = AIMSoln[rssmodel,{},simdegree,AIMZeroTol,AIMPrecision][[Flatten[Map[Position[AIMVarNames[rssmodel], #] &,
													simvars]]]] ;

soln = Chop[N[soln /.Last[AIMGenericArgs[rssmodel]]->1 /. momSubs]] /.
	Thread[Drop[AIMGenericArgs[rssmodel],-Length[AIMShocks[rssmodel]]-1] -> Map[inputvar[[#]]&,lagvarpos]] /.
	Thread[Take[AIMGenericArgs[rssmodel],{-Length[AIMShocks[rssmodel]]-1,-2}] -> 
							Map[inputshock[[#]] &, Range[Length[AIMShocks[rssmodel]]]]] ;

iterator = Apply[Function,{ {inputvar,inputshock}, soln} ] ;


(********This is probably what we want!!!!!!!!!!!!!!!!****************)

syndata = Transpose[FoldList[iterator, Table[0,{Length[simvars]}], shocks]] ;

synpricebond = Flatten[syndata[[Flatten[Position[simvars, pricebond]]]]] + (pricebondAIMSS /.ss) ;
synytm = Flatten[syndata[[Flatten[Position[simvars, ytm]]]]] + (ytmAIMSS /.ss) ;
synytmrn = Flatten[syndata[[Flatten[Position[simvars, ytmrn]]]]] + (ytmAIMSS /.ss) ;
synInt = Flatten[syndata[[Last[Position[simvars,Int]]]]] + (IntAIMSS /.ss) ;
syntp = Flatten[syndata[[Flatten[Position[simvars, termprem]]]]] + (termpremAIMSS /.ss) ;

meantp = Mean[syntp] ; (* mean term prem in annualized basis points *)

stddevs = Map[StandardDeviation,syndata[[reportvarpos]]] * {100, 100, 100, 100, 400, 400, 400, 1, 1}
		(* real variable std devs in percentage points, not logs; inflation std devs in
			annualized percentage points; interest rate std devs in annualized basis points *)

meanslope = (Mean[synytm] - 400*Mean[synInt]) * 100 ; (* slope in annualized basis points *)
meanslopern = (Mean[synytmrn] - 400*Mean[synInt]) * 100 ; (* slope in annualized basis points *)
stdslope = StandardDeviation[synytm - 400*synInt] * 100 ;

synehpr = (((consoldelta /.parametervalues/.ss) *Drop[Exp[synpricebond],1] + Drop[Exp[synInt],-1] *1/100) /
		Drop[Exp[synpricebond],-1] - Drop[Exp[synInt],-1]) * 40000 ; (* ehpr in annualized basis points *)
meanehpr = Mean[synehpr] ;
stdehpr = StandardDeviation[synehpr] ;

synDpricebond = Drop[synpricebond,1] - Drop[synpricebond,-1] ;
cscoeff = Covariance[-synDpricebond,Drop[synytm/400 - synInt,-1]] / Variance[synytm/400 - synInt] ;



Print["simulated mean tp = ", meantp]
Print["simulated mean slope = ", meanslope]
Print["simulated mean ehpr = ", meanehpr]
Print["simulated std dev tp = ", stdtp]
Print["simulated std dev slope = ", stdslope]
Print["simulated std dev ehpr = ", stdehpr]

Export["testpapermodel.txt",syndata,"Table"]



(*
  delta -> .02,
  phi -> 2,
  habitsize -> .66,
  chi0 -> (* 4.73921522897550247856008201926876616780052304558537869572838202224763 *)
    Exp[wrealAIMSS +MUcAIMSS], (* normalize L in the nonstochastic steady state to be 1; note the use of *AIMSS
				variables here to accomplish this.  The AIMSS suffix denotes the nonstochastic
				steady-state value of a variable, as defined by the AIMSS[_] function call. *)
  chi -> 1.5,

  psi -> 0 *YBar, (* quadratic adjustment costs to labor *)
  theta -> .2,
  xi -> .75,

  ABar -> 1,
  GBar -> .17 *YBar,
  KBar -> 10 *YBar,
  IBar -> delta *KBar,
  YBar -> Exp[YAIMSS], (* nonstochastic steady-state level of Y *)
  piBar -> 0,

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

  rhoa -> .9,
  rhog -> .9,
  rhopistar -> 0,
  rhoinflavg -> .7,
*)