3. Hit the ground running: an introductory example

This chapter mirrors the canonical introductory example used in the legacy documentation. The same model.rs file and the same calibration / priors / sampling code are used; only the constructor (dsge_model(...) instead of rise(...)) and the named-argument form of set / estimate differ. Reading this chapter against the legacy Teaser Example makes the stable-to-modern syntax delta concrete in one place.

3.1. The model economy 

A three-equation New Keynesian model. The first equation is an Euler-type IS curve relating output to the real interest rate.

(1)\[ \beta E_{t}\frac{r_{t}}{\pi _{t+1}}\left[ \frac{1-\frac{\kappa }{2}\left( \pi _{t}-\pi _{t-1}^{\alpha }\pi ^{1-\alpha }\right) ^{2}}{1-\frac{\kappa }{2} \left( \pi _{t+1}-\pi _{t}^{\alpha }\pi ^{1-\alpha }\right) ^{2}}\right] \frac{Y_{t}}{Y_{t+1}}=1\]

The second equation is a nonlinear Phillips curve with price adjustment costs.

(2)\[\begin{split} \begin{array}{ccc} 0 & = & \chi _{0}\eta _{t}\left[ 1-\frac{\kappa }{2}\left( \pi _{t}-\pi _{t-1}^{\alpha }\pi ^{1-\alpha }\right) ^{2}\right] \left( \frac{Y_{t}}{Z_{t}}\right) ^{1+\chi } \\ & & +\left( 1-\eta _{t}\right) -\kappa \left( \pi _{t}-\pi _{t-1}^{\alpha }\pi ^{1-\alpha }\right) \pi _{t} \\ & & +\beta \kappa E_{t}\left( \pi _{t+1}-\pi _{t}^{\alpha }\pi ^{1-\alpha }\right) \pi _{t+1}\left[ \frac{1-\frac{\kappa }{2}\left( \pi _{t}-\pi _{t-1}^{\alpha }\pi ^{1-\alpha }\right) ^{2}}{1-\frac{\kappa }{2}\left( \pi _{t+1}-\pi _{t}^{\alpha }\pi ^{1-\alpha }\right) ^{2}}\right] \end{array}\end{split}\]

The third equation is a monetary policy reaction function with a state-dependent inflation reaction.

(3)\[ \frac{r_{t}}{r}=\left( \frac{r_{t-1}}{r}\right) ^{\rho }\left[ \left( \frac{ \pi _{t}}{\pi _{\ast }}\right) ^{\psi _{\pi }\left( s_{t}\right) }\right] ^{\left( 1-\rho \right) }\exp \left( \sigma _{r}\varepsilon _{r,t}\right)\]

The fourth equation is the cost-push shock process.

(4)\[ \log \eta _{t}=\left( 1-\rho _{\eta }\right) \log \eta +\rho _{\eta }\log \eta _{t-1}+\sigma _{\eta }\varepsilon _{\eta ,t}\]

The fifth equation is the technology shock process.

(5)\[ \log \left( \frac{Z_{t}}{Z_{t-1}}\right) =\left( 1-\rho _{z}\right) \mu _{z}\left( s_{t}\right) +\rho _{z}\log \left( \frac{Z_{t-1}}{Z_{t-2}}\right) +\sigma _{z}\varepsilon _{z,t}\]

Regime switching enters through \(\psi_{\pi}\) (the central bank’s inflation reaction) and \(\mu_{z}\) (the growth rate of technology), switching in lockstep on a chain we name snk. Uppercase variables are non-stationary; lowercase variables are stationary.

3.2. The model file 

The same model.rs file used in the legacy documentation. The file lives next to this chapter; full listing:

%% Simple New Keynesian model: nonstationary

@endogenous(log)	Y "Output", R "interest rate", PAI "Inflation"
Z "Technology", ETA "Cost-push shock process", DY "Output growth"

@exogenous EZ "Technology shock", ER "Monetary policy shock", EETA "Cost-push shock"

@parameters	beta "discount factor", kappa "Price adjust. cost coef",
alpha "indexation on past inflation", eta "steady-state elasticity of subst. across goods",
chi0 "scale of disutility of labor", chi "Inverse frisch elasticity of subst.",
rho "interest rate smoothing", rhoz "persistence: technology", rhoeta "persistence: cost-push",
sigr "scale: mon. pol. shock", sigz "scale: technology shock",
sigeta "scale: cost-push shock", paitarg "Inflation target"
snk_tp_1_2, snk_tp_2_1

@parameters(snk,2) muz "growth-rate of technology", psi_pai "reaction function: inflation"

@observables PAI, R, DY

@model

	"IS curve"
	beta*R{t}/PAI{t+1}*[(1-kappa/2*(PAI{t}-PAI{t-1}^alpha*PAI{stst}^(1-alpha))^2)/...
		(1-kappa/2*(PAI{t+1}-PAI{t}^alpha*PAI{stst}^(1-alpha))^2)]*Y{t}/Y{t+1} = 1;

	"Phillips curve"
	0 = chi0*ETA*[1-kappa/2*(PAI{t}-PAI{t-1}^alpha*PAI{stst}^(1-alpha))^2]*(Y/Z)^(1+chi)+...
		(1-ETA{t})-kappa*(PAI{t}-PAI{t-1}^alpha*PAI{stst}^(1-alpha))*PAI{t}+...
		beta*kappa*(PAI{t+1}-PAI{t}^alpha*PAI{stst}^(1-alpha))*PAI{t+1}*...
		[(1-kappa/2*(PAI{t}-PAI{t-1}^alpha*PAI{stst}^(1-alpha))^2)/...
		(1-kappa/2*(PAI{t+1}-PAI{t}^alpha*PAI{stst}^(1-alpha))^2)];

	"Monetary policy reaction function"
	R{t}/R{stst} = (R{t-1}/R{stst})^rho*[(PAI{t}/paitarg)^psi_pai]^(1-rho)*exp(sigr*ER{t})
	# PAI{t}=paitarg;

	"Cost-push shock process"
	log(ETA{t}) = (1-rhoeta)*log(eta) + rhoeta*log(ETA{t-1}) + sigeta*EETA{t};

	"Technology shock process"
	log(Z{t}/Z{t-1}) = (1-rhoz)*log(muz) + rhoz*log(Z{t-1}/Z{t-2}) + sigz*EZ{t};

	"Output growth"
	DY = Y{t}/Y{t-1};
	

The two distinguishing declarations:

@endogenous(log) Y "Output", R "interest rate", PAI "Inflation", Z "Technology", ETA "Cost-push shock process", DY "Output growth" – every endogenous variable is log-linearised.
@parameters(snk,2) muz "growth-rate of technology", psi_pai "reaction function: inflation" – the two parameters that switch, on a two-state chain named snk.

The full file is model.rs.

3.3. Building the model object 

m = dsge_model('model');

(The legacy documentation used m = rise('model') or m = dsge('model'). The modern toolbox replaces the class-named constructor with shape-specific factories; dsge_model is the DSGE factory.)

Without further options, RISE computes and stores first-order symbolic derivatives suitable for a first-order perturbation.

3.4. Parameterizing the model 

p                = struct();
p.beta           = 0.99;
p.kappa          = 161;
p.paitarg        = 1.02^0.25;
p.alpha          = 0.5;
p.eta            = 6;
p.chi            = 0.7;
p.chi0           = 1;
p.psi_pai_snk_1  = 2.5;
p.psi_pai_snk_2  = 0.9;
p.rho            = 0.7;
p.rhoz           = 0.75;
p.rhoeta         = 0.75;
p.muz_snk_1      = 1.04^.25;
p.muz_snk_2      = 1.01^.25;
p.sigz           = 0.05;
p.sigeta         = 0.05;
p.sigr           = 0.05;
p.snk_tp_1_2     = 1 - 0.95;
p.snk_tp_2_1     = 1 - 0.9;

m = set(m, parameters = p);

Switching parameters carry the chain name and state in their name: psi_pai_snk_1 is \(\psi_{\pi}\) in state 1 of chain snk; snk_tp_1_2 is the transition probability from state 1 to state 2 of chain snk.

3.5. Solving the non-stationary model 

The model is non-stationary; we let RISE solve it on the balanced-growth path with the Maih-Waggoner perturbation strategy and the matching Newton solver:

m = solve(m, ...
    solve_perturbation_type = 'mw', ...
    solve_bgp               = true, ...
    solver                  = 'mn');

print_solution(m);

3.6. A one-shot `stoch_simul` run 

info = stoch_simul(m);

info carries simulated data, moments, IRFs, autocorrelations, skewness and kurtosis.

3.7. Impulse responses 

myirfs = irf(m);
quick_irfs(m, myirfs, get(m, 'endo_list(original)'));

The output is a struct of ts keyed first by shock name then by endogenous variable; endo_list(original) excludes the auxiliary variables RISE introduces internally.

3.8. 5th-order perturbation 

m5 = solve(m, ...
    solve_order            = 5, ...
    solve_derivatives_type = 'automatic');

print_solution(m5);

The model was originally parsed for first-order symbolic derivatives; for fifth order we switch the derivative engine to automatic (algorithmic) differentiation. Alternatively, re-parse the model with max_deriv_order = 5 to get symbolic derivatives up to order 5.

3.9. Data from FRED 

xrange = '1960Q1:2022Q3';
rawdb  = fetch_fred({'BPCCRO1Q156NBEA', ...
                     'BOGZ1FL072052006Q', ...
                     'NAEXKP01USQ657S'});

db        = struct();
db.PAI    = 1 + rawdb(1).series(xrange)/100;
db.R      = 1 + rawdb(2).series(xrange)/100;
db.GROWTH = 1 + rawdb(3).series(xrange)/100;

These are the canonical PCE-inflation, federal funds rate, and US-output series used in the documentation example.

3.10. Priors 

The legacy chapter sets priors in quantile form. We do the same. Several priors illustrate the four-, five-, and six-slot variants (see Estimation for the full parametrisation taxonomy):

priors = struct();

% 4-slot quantile form
priors.kappa = {p.kappa, 5,    20,   'gamma(.9)'};
priors.alpha = {p.alpha, 0.05, 0.948,'beta(.9)'};
priors.chi   = {p.chi,   1.5,  3,    'gamma(.9)'};

% 6-slot generalized-beta form (extra hyperparameters)
priors.eta       = {p.eta,        3,        8,        1,    12,        'beta(.9)'};
priors.muz_snk_1 = {p.muz_snk_1,  1.02^.25, 1.05^.25, 1,    1.07^.25,  'beta(.9)'};
priors.muz_snk_2 = {p.muz_snk_2,  1.01^.25, 1.03^.25, 1,    1.04^.25,  'beta(.9)'};

% 4-slot quantile form
priors.rho    = {p.rho,    0.05, 0.948, 'beta(.9)'};
priors.rhoz   = {p.rhoz,   0.05, 0.948, 'beta(.9)'};
priors.rhoeta = {p.rhoeta, 0.05, 0.948, 'beta(.9)'};

% 4-slot + hard truncation (6-slot)
priors.sigz   = {p.sigz,   0.0005, 1.0, 'sichisq(.9)', 0, 3};
priors.sigeta = {p.sigeta, 0.0005, 1.0, 'sichisq(.9)', 0, 3};
priors.sigr   = {p.sigr,   0.0005, 1.0, 'sichisq(.9)', 0, 3};

% 4-slot, transition probabilities
priors.snk_tp_1_2 = {p.snk_tp_1_2, 0.01, 0.411, 'beta(.9)'};
priors.snk_tp_2_1 = {p.snk_tp_2_1, 0.01, 0.411, 'beta(.9)'};

% 6-slot generalized-beta, switching reaction
priors.psi_pai_snk_1 = {p.psi_pai_snk_1, 1.1, 2.5, 1,   3,   'beta(.9)'};
priors.psi_pai_snk_2 = {p.psi_pai_snk_2, 0.5, 1.1, 0.1, 1.5, 'beta(.9)'};

3.11. Visualizing the priors 

plotOpts             = struct();
plotOpts.prior_trunc = 2.6e-3;

rdist.plot(priors, plotOpts, [], struct('linewidth', 2));

3.12. Posterior maximisation 

mest = estimate(m, ...
    data             = db, ...
    estim_priors     = priors, ...
    kf_init_variance = 1, ...
    kf_presample     = 10, ...
    optimizer        = 'bee_gate');

kf_init_variance = 1 – initial variance of every endogenous variable, set to 1 because the model is non-stationary and does not have a finite unconditional variance.
kf_presample = 10 – the first 10 observations are not included in the likelihood; this is a defensive complement to kf_init_variance.
optimizer = 'bee_gate' – the artificial bee colony global-search optimizer, robust on this calibration.

3.13. Posterior simulation 

We use the random-walk Metropolis-Hastings sampler from the modern +rsamplers package. pull_objective returns a value to minimize; rsamplers.rwmh takes a value to maximize, so we negate the objective before passing it in:

[objective, lb, ub, x0, SIG] = pull_objective(mest);

scale           = 0.0794;
myOpts          = struct();
myOpts.tunedCov = scale * SIG;
myOpts.N        = 20000;

energy  = @(varargin) -objective(varargin{:});
results = sample(rsamplers.rwmh(energy, x0, lb, ub, myOpts));

results{1}.stats

3.14. Marginal data density via the bridge estimator 

RISE ships nine ways of computing the marginal data density. For the Meng-Wong (1996) bridge estimator:

mddobj = mdd(results, energy, lb, ub, [], [], true);
bridge(mddobj, true, mdd.global_options);

3.15. The full driver 

The complete driver script is model_driver.m shipped alongside this chapter. Running it end to end requires network access for the FRED fetch.

The legacy documentation ships the same model file and a driver with the same content under the legacy syntax. Diffing the two drivers side by side is the cleanest way to see the stable-to-modern syntax delta in practice.