simulated moments

[otb]

When the moments are calculated from simulations, they are calculated from a single simulation (of size=periods), correct? Is there quick way to take the average of multiple simulations to get the moments? Thanks.

[otb]

Let me ask some related questions, and maybe I can get answers to my previous questions as well.

Should the moments generated from a single simulation approach the true theoretical moments as we increase periods?

If so, why do some papers generate moments from multiple simulations each of which have a small number of periods (typically equal to the number of periods in the data)?

I guess this is related to the distinction between consistency vs. unbiasedness of an estimator. So the two methods (single simulation with infinte periods vs. the mean of infinite simulations with finite periods) should give the same answer, right?

Thanks for any help.

[StephaneAdjemian]

Hi otb,

Let me ask some related questions, and maybe I can get answers to my previous questions as well.

Should the moments generated from a single simulation approach the true theoretical moments as we increase periods?

Yes, this should be true if the DGP is stationary.

If so, why do some papers generate moments from multiple simulations each of which have a small number of periods (typically equal to the number of periods in the data)?

Let us consider a simple example. Suppose that our DGP is the following AR(1) stochastic process:

Code: Select all

y_t = c + \rho y_{t-1} + \epsilon_t


where |\rho|<1 and {\epsilon_t} is a mean zero white noise with variance \sigma^2, so that {y_t} is an asymptotically second order stationary. One can easily show that the asymptotic first and second order moments (the theoretical moments in Dynare words) are:

Code: Select all


   E_{\infty}[y_t] = \frac{c}{1-\rho}

   V_{\infty}[y_t] = \frac{\sigma^2}{1-\rho^2}



Now if we suppose that there exists an initial condition, say, y_0 = 0, the moments will be different from the asymptotic moments because they depend on the initial condition. To understand the reason we just have to note that :

Code: Select all


   E[y_1] = c + \rho E[y_0] = (1 + \rho) \times c

   V[y_1] = \rho^2 V[y_0] + \sigma^2 = (1+\rho^2) \times \sigma^2



Code: Select all


   E[y_2] = c + \rho E[y_1] = (1 + \rho + \rho^2) \times c

   V[y_2] = \rho^2 V[y_1] + \sigma^2 = (1+\rho^2+\rho^4) \times \sigma^2



and more generally

Code: Select all


   E[y_t] =  (1 + \rho + ... +\rho^t) \times c

   V[y_t] =  (1+\rho^2+ ... + \rho^{2t}) \times \sigma^2



or equivalently:

Code: Select all


   E[y_t] =  \frac{1-\rho^{t-1}}{1-\rho} \times c

   V[y_t] =  \frac{1-\rho^{2(t-1)}}{1-\rho^2} \times \sigma^2



Obviously when t goes to infinity these moments converge to the asymptotic moments. The speed of convergence depends on the value of the autoregressive parameter. The closer to one is \rho, the lower will be the convergence to the asymptotic moments.

One can check that the formulas for the expectation will be affected if we change the initial condition. For instance, if we set y_0 equal to E_{\infty}[y_t] instead of 0 we can see that the expectation becomes time invariant. More generally if the initial condition is a random variable (with non zero variance) the formulas for the second order moment will also be affected.

To sum-up, simulated moments have to be different from asymptotic moments because of the influence of the initial condition on the moments (and also, but that is obvious, because of sampling issues).

Averaging over different path simulations starting from different initial conditions is a way to overcome this difficulty. Ideally we should randomly select initial conditions in a distribution with expectation E_{\infty}[y_t] and
variance V_{\infty}[y_t], but this is not possible since we don't know these moments in the first place.

A single simulation with an infinite number of periods will give results different from the mean of an infinite number of simulations of length t<\infty (because the process is only asymptotically second order stationary) except if the initial conditions are obtained from the ergodic distribution of the stochastic process (which is unknown). A single simulation with an infinite number of periods will converge to E_{\infty}[y_t]
and V_{\infty}[y_t] in probability, while the mean of an infinite number of simulations of length t will converge in probability to E[y_t] and V[y_t] (these moments will depend on the assumption made about the distribution of the initial condition).


Best Regards,
Stéphane.

[otb]

Thanks for the detailed response Stephane. This was very helpful. Since I want asymptotic moments, I should be fine.

But in any case, does Dynare conduct multiple simulations? I don't think so, but I just want to make sure and also see if somebody had a quick fix to actually make Dynare do that (perhaps write a loop in the mod file).
Thanks again.

otb

[StephaneAdjemian]

Thanks for the detailed response Stephane. This was very helpful. Since I want asymptotic moments, I should be fine.

But in any case, does Dynare conduct multiple simulations?

No. But you can code a matlab loop around the stoch_simul command.

Best, Stéphane.

[quiko]

Hi,
I have a related question. I wrote a loop around stoch_simul(hp_filter = 1600, order = 1) over a continuum of values for the parameter of my interest. I am interested in the theoretical moments of each parameter's iteration. They are printed in the command window, BUT I would like to have them stored somewhere for each iteration: those ones contained in oo_ refer only to the last iteration, but I would like to have them stored for each iteration. How could I solve this? Thanks for your help.

This is part of the loop I have written:
[...]
steady;
gammas = 0.3:0.05:0.5; // weights of risk
for s = 1:length(gammas);
gamma_y = gammas(s);
stoch_simul(hp_filter = 1600, order = 1, nograph);
if info;
disp(['Computation fails for gamma_y = ' num2str(gamma_y)]);
end;
// std deviation here

[jpfeifer]

Add an else clause where you save the contents of oo_, i.e. something along the lines of

Code: Select all

if info;
disp(['Computation fails for gamma_y = ' num2str(gamma_y)]);
else
moments{s,1}=oo_;
end;

[quiko]

Thanks for the reply, it works. One last question: If I want to extrapolate from this structure only the standard deviation of each iteration, how the code would look like (converting into .mat) ? Thanks.

[jpfeifer]

It depends what exactly you need. Variances are stored in oo_.variance, so

Code: Select all

if info;
disp(['Computation fails for gamma_y = ' num2str(gamma_y)]);
else
moments(:,s)=diag(oo_.var);
end;

[quiko]

I get "Reference to non-existent field 'variance', but I changed the command to moments{s,:} = oo_.var; which gives me as a result a <5X1 cell>. Each cell is a <1X1 struct> inside which there is the variance of each variable. Ideally, what I exactly need is a way to put these variances for each into a matrix M(n,s).

Hope I have been clear. Thank you again for your useful help.

[jpfeifer]

My mistake. I updated the code above. Note the use of round brackets instead of curly ones.

[quiko]

Thanks a lot, it's perfect.

[quiko]

It depends what exactly you need. Variances are stored in oo_.variance, so

Code: Select all

if info;
disp(['Computation fails for gamma_y = ' num2str(gamma_y)]);
else
moments(:,s)=diag(oo_.var);
end;

Hi, what if instead of the variance I want to get the covariance matrix? I tried with the following command:

covars(:,s)=triu(oo_.var);

but I get an error message due to dimension mismatch.

Thank you very much.

[jpfeifer]

triu provides a matrix, but the left-hand-side expects a vector. That cannot work. You want something like

Code: Select all

covars=triu(oo_.var);

[quiko]

I think I have not been clear enough.
Recalling that n is the number of variables and s is the number of iterations, in my first question I wanted to see how the variance changes for each iteration. So, the command: moments(:,s)=diag(oo._var) did the right job by providing as output a nXs matrix. What I exactly need now is to have the covariance, but I want it for each variable with respect to all the others AND for each iteration. The command: covars=triu(oo._var) does't do the job because I get a nX1 vector.

Maybe what I need is the first raw of (oo._var), which should be the covariance of the first variable with respect to all the others. Then I could repeat manually the second raw, the third, and so on. But I don't know how to get this first raw from (oo_.var) instead of its diagonal, how to get something like: covars(:,s)=first_raw(oo._var).

Would it be correct the following command?: covars(:,s)=(oo_.var(1,:))

Thank you for your replies.