Impulse responses with multiple shocks

I have a simple RBC model with time-varying risk of the form

// stochastic process
sigma_e = rho_sigma_e*sigma_e(-1) + (1-rho_sigma_e)*mean_sigma_e + sigma_eta*eta;
a = rho_a*a(-1) + sigma_e*e;

where eta and e are exogenous random variables and the neoclassical production function is

y = exp(a) * k^alpha * l^(1-alpha)

I would like to use Dynare to generate impulse response functions to risk shocks (i.e. impulse responses to eta). I am performing a third-order approximation of the model so the policy functions depend on sigma_e. Thus a shock to risk (sigma_e) should generate temporary divergence from steady state, even if all draws of the exogenous variable e are zero. This is precisely what I would like to do. However, this is not what Dynare does for the impulse responses to eta. Instead, it seems to draw values of e from its given distribution. Here is what the IRF to eta looks like:

You can see that TFP (a) is receiving non-zero shocks. What I want is for TFP to remain at its steady state value (i.e. the shock e should be zero throughout the impulse response) so that the responses of the endogenous variables to risk shocks are clear. Is there a way to do this? I have attached my complete mod file for reference. Thanks in advance.
one_country_RBC_time_var_vol.mod (2.07 KB)

I have a similar problem. When I shock the volatility (sigma_e), tfp (a) moves around (I assume because shocks e are also being drawn to simulate the IRFs). However, when I shock tfp directly (an e shock), there is no reaction whatsoever in volatility (sigma_e response is exactly zero). Which makes me feel that shocks eta are not being drawn in this case. Why is does this asymmetry happen? I checked the Dynare manual, but the explanation on Monte-Carlo simulations to compute IRFs for higher order approximations is very generic.

Any help is very much appreciated.

Have you figured out how to do this yet, jsteinberg?

No I have not. I would still like to know how to do it. Alternatively, could someone explain to me how to manually run simulations using the third-order decision rules? It’s easy for 1st order, but it’s not so clear for higher order approximations.

I figured out how to do this by manually calling the simult_ function. One of its inputs is the sequence of exogenous shocks, so I just constructed that sequence like I wanted and called the function. The output is the time series of all the endogenous variables after impulse.

Can you show me how you did that?


No news about it? I am doing volatility shocks as well but I am facing the same problem you posted.
Any help would be useful.

In Dynare, we compute the IRF for shock e, for order=2,3 in the following manner:

  1. draw a series of random shocks for 140 periods
  2. Perform simulation Y1
  3. Add one standard deviation to the simulated series for shock e in period 101
  4. Perform simuluation Y2
  5. The IRF for this experiment is equal to Y2-Y1
  6. Results obtained in 5) is affected by idiosyncratic shocks other than the one impulse in period 101. In order to average over the effect of these idiosyncratic shocks, we perform 50 replications of steps 1) to 5) and report the average.

-Initial 100 warming periods is the default value of option DROP
-40 periods for the IRF is the default value of option IRF
-50 replications is the default value of option REPLIC

By following this procedure, we average over the influence of the state point at which the impulse occurs and of future shocks in periods after the impulse. So, we report an average IRF

In the first example in this thread, an impulse on eta has an effect on sigma_e and therefore on a and the other variables of the model. It also creates numerical noise that is very difficult to eliminate by increasing the number of replications.

When jsteinberg says

I think that what you want is the effects on a shock on eta, conditional on the realization of e being always zero. We don’t have tool to do that, but, as mentioned by jsteinberg in the thread, it can be done by using function simult_.

I attach an exemple on how to show the IRF for consumption


one_country_RBC_time_var_vol_MJ.mod (2.2 KB)

Thank you very much

Thanks for good response with multiple stocks.

This was very useful explanation!

I would, however, like to ask a further clarification.
Using the example above *_MJ.mod I tried to do the same exercise at order=3.
The IRFs somewhat differ from what we get at order=2. I attach the figures. In particular what confuses me is:
1.) at order=2, the shock to volatility increases volatility then it returns to its steady state value. However, at order=2, the volatility shock should not have the impact on the path of aggregate variables (only the constant in the policy function).
2.) At order=3, I expect the same behaviour of volatility as at order=2. However, it seems as it drops (after a positive shock) and remains at lower level. This obviously results in strange path of consumption.

The results must come from the line change of simult_ command. With the last argument being order as far as I understand (function y_=simult_(y0,dr,ex_,iorder)].

I changed
irf2 = simult_(oo_.steady_state,oo_.dr,ex_,2);
irf3 = simult_(oo_.steady_state,oo_.dr,ex_,3);

Probably I misinterpreted something. Any ideas or suggestions on this point?

Modified *_MJ.mod code:

var a, sigma_e, c, l, x, k, y;
varexo e, eta;

parameters alfa, betta, gamma, mu, delta, phi, rho_a, rho_sigma_e, mean_sigma_e, sigma_eta;

alfa = 0.36; // capital’s share of output
betta = 0.99; // discount factor
gamma = -1; // RRA = 1 - gamma
mu = 0.32; // exponent on leisure in Cobb-Douglas part of utility
delta = 0.025; // depreciation
phi = 3; // capital adjustment factor
rho_a = 0.95; // persistence of TFP
rho_sigma_e = 0.97; // persistence of std dev of TFP innovations
mean_sigma_e = 0.008; // unconditional mean of std dev of TFP innovations
sigma_eta = 0.0003; // std dev of innovations to TFP volatility


// output
y = exp(a)(k(-1)^alfa)(l^(1-alfa));

// labor-leisure

// budget constraint
c = y - x;

// firm FOC for investment
(((c^mu)((1-l)^(1-mu)))^(gamma-1))(mu*(c^(mu-1))((1-l)^(1-mu)))=(1-2phi*(x/k(-1)-delta))betta(((c(+1)^mu)((1-l(+1))^(1-mu)))^(gamma-1))(mu*(c(+1)^(mu-1))((1-l(+1))^(1-mu)))( alfay(+1)/k + (1-delta + phi((x(+1)/k)^2 - delta^2))/(1-2phi(x(+1)/k-delta)) );

// law of motion for capital
k = (1-delta)k(-1) + x - phik(-1)*(x/k(-1)-delta)^2;

// stochastic process
sigma_e = rho_sigma_esigma_e(-1) + (1-rho_sigma_e)mean_sigma_e + sigma_etaeta;
a = rho_a
a(-1) + sigma_e*e;



sigma_e = mean_sigma_e;
x = delta
c = y - x;



var e = 1;
var eta = 1;

stoch_simul(order=3, periods=10000, replic=100);

ex_ = zeros(10,2);
ex_(1,2) = 1;

ex_(1,2) = 1;
irf2 = simult_(oo_.steady_state,oo_.dr,ex_,2);
irf3 = simult_(oo_.steady_state,oo_.dr,ex_,3);




close all;

legend(‘A 2nd ord’)
legend(‘A 3rd ord’)

legend(‘sigmae 2nd ord’)
legend(‘sigmae 3rd ord’)

legend(‘cons 2nd ord’)
legend(‘cons 3rd ord’)


The solution to my question is in the post: Using simult_ with third order approximation
Then, everything seems ok.
Please, correct me if you spot that it is not…