(* Equations of the model *)
eqns={
Y[t] == A[t] *K[t-1]^alpha,
Log[A[t]] == rho *Log[A[t-1]] + eps[a][t],
K[t] == (1-delta) *K[t-1] + Inv[t], (* Mathematica treats "I" as Sqrt[-1] *)
Y[t] == C[t] + Inv[t],
C[t]^-gamma == beta *(1+r[t+1]) *C[t+1]^-gamma,
r[t] == alpha *A[t] *K[t-1]^(alpha-1) - delta,
Welf[t] == C[t]^(1-gamma) /(1-gamma) + beta *Welf[t+1]
}
(* variables to transform to logs and then approximate *)
logvars = {A, C, K, Y}
logrules = Map[#[x_]->E^(#[x])&, logvars]
(* parameter values *)
parametervals={
alpha->0.3`70, (* 70 digits of precision: in general, you need to *)
beta->0.99`70, (* input more precision than you desire for the output. *)
gamma->1.1`70, (* Also, the larger the model and the higher the order *)
delta->0.1`70, (* of approximation, the more precision you must input *)
rho->0.8`70 (* to achieve a given desired precision of output. *)
}
(* substitute log transformation rules and parameter values into equations: *)
sgmodel = eqns /.logrules //.parametervals
(* complete variable list (Mathematica places in alphabetical order):
{A, C, Inv, K, r, Welf, Y}
*)
(* Now find the steady state to about 50 significant digits. *)
(* To truly get 50 significant digits of precision, one should set the
AIMZeroTol parameter to 10^-50. I have not done this here because I
want to avoid giving users the misimpression that this is what they
should do for their models in general. As discussed in the comments to
perturbationAIMsimple.m, the AIMZeroTol parameter should be thought of
essentially as an "economic zero" rather than a "numerical zero" (the
latter is left to the internals of Mathematica to handle appropriately).
Since we're using arbitrary-precision numbers here, we have more freedom
to set AIMZeroTol conservatively, but values much smaller than 10^-20 are
probably not very meaningful economically.
*)
AIMZeroTol=10^-20 ;
(* The AIMPrecision option essentially determines the precision of the
output of AIMSS; it must be less than the precision of the inputs or you
will get a warning. You must calculate the steady state to greater
precision than you desire for the output of AIMSeries.
*)
AIMSS[sgmodel, AIMPrecision->65] ;
(* Find first- and second-order coefficients to 50 significant digits *)
Print[AIMSeries[sgmodel,1] //TableForm]
AIMSeries[sgmodel,2] //TableForm
(* Truncate the number of digits reported in the answer, if you prefer
(the following reports only the first 6 significant digits, even though
the answer was computed to about 50):
NumberForm[TableForm[AIMSeries[sgmodel,2]],6]
*)
(* Check the accuracy of the output (reports number of significant digits
to the right of decimal point in the answer):
Accuracy[AIMSeries[sgmodel,2]]
*)