# Simulated moments

**URL:** <https://forum.dynare.org/t/simulated-moments/670>\
**Category:** Dynare help (legacy posts)\
**Created:** [11 May 2008 01:36 UTC](https://forum.dynare.org/t/simulated-moments/670 "2008-05-11T01:36:23Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![otb](https://forum.dynare.org/letter_avatar_proxy/v4/letter/o/46a35a/32.png) [@otb](https://forum.dynare.org/u/otb)\
**Post date:** [11 May 2008 01:36 UTC](https://forum.dynare.org/t/simulated-moments/670/1 "2008-05-11T01:36:23Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![otb](https://forum.dynare.org/letter_avatar_proxy/v4/letter/o/46a35a/32.png) [@otb](https://forum.dynare.org/u/otb)\
**Post date:** [13 May 2008 02:43 UTC](https://forum.dynare.org/t/simulated-moments/670/2 "2008-05-13T02:43:35Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stepan-a](https://forum.dynare.org/user_avatar/forum.dynare.org/stepan-a/32/6652_2.png) [@stepan-a](https://forum.dynare.org/u/stepan-a)\
**Post date:** [13 May 2008 10:26 UTC](https://forum.dynare.org/t/simulated-moments/670/3 "2008-05-13T10:26:57Z")

</div>

Hi otb,

[quote=“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? [/quote]

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

> [@otb](#):
>
> 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:

```auto
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:

```auto

   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 :

```auto

   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

```

```auto

   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

```auto

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

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

```

or equivalently:

```auto

   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.

---

<div class="post-metadata">

**Author:** ![otb](https://forum.dynare.org/letter_avatar_proxy/v4/letter/o/46a35a/32.png) [@otb](https://forum.dynare.org/u/otb)\
**Post date:** [13 May 2008 13:39 UTC](https://forum.dynare.org/t/simulated-moments/670/4 "2008-05-13T13:39:37Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![stepan-a](https://forum.dynare.org/user_avatar/forum.dynare.org/stepan-a/32/6652_2.png) [@stepan-a](https://forum.dynare.org/u/stepan-a)\
**Post date:** [15 May 2008 10:21 UTC](https://forum.dynare.org/t/simulated-moments/670/5 "2008-05-15T10:21:43Z")

</div>

[quote=“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? [/quote]

No. But you can code a matlab loop around the stoch\_simul command.

Best, Stéphane.

---

<div class="post-metadata">

**Author:** ![quiko](https://forum.dynare.org/letter_avatar_proxy/v4/letter/q/49beb7/32.png) [@quiko](https://forum.dynare.org/u/quiko)\
**Post date:** [10 April 2014 22:34 UTC](https://forum.dynare.org/t/simulated-moments/670/6 "2014-04-10T22:34:16Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jpfeifer](https://forum.dynare.org/user_avatar/forum.dynare.org/jpfeifer/32/5044_2.png) [@jpfeifer](https://forum.dynare.org/u/jpfeifer)\
**Post date:** [11 April 2014 05:25 UTC](https://forum.dynare.org/t/simulated-moments/670/7 "2014-04-11T05:25:42Z")

</div>

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

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

---

<div class="post-metadata">

**Author:** ![quiko](https://forum.dynare.org/letter_avatar_proxy/v4/letter/q/49beb7/32.png) [@quiko](https://forum.dynare.org/u/quiko)\
**Post date:** [11 April 2014 09:53 UTC](https://forum.dynare.org/t/simulated-moments/670/8 "2014-04-11T09:53:51Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![jpfeifer](https://forum.dynare.org/user_avatar/forum.dynare.org/jpfeifer/32/5044_2.png) [@jpfeifer](https://forum.dynare.org/u/jpfeifer)\
**Post date:** [11 April 2014 11:48 UTC](https://forum.dynare.org/t/simulated-moments/670/9 "2014-04-11T11:48:18Z")

</div>

It depends what exactly you need. Variances are stored in oo\_.variance, so

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

---

<div class="post-metadata">

**Author:** ![quiko](https://forum.dynare.org/letter_avatar_proxy/v4/letter/q/49beb7/32.png) [@quiko](https://forum.dynare.org/u/quiko)\
**Post date:** [11 April 2014 13:03 UTC](https://forum.dynare.org/t/simulated-moments/670/10 "2014-04-11T13:03:16Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![jpfeifer](https://forum.dynare.org/user_avatar/forum.dynare.org/jpfeifer/32/5044_2.png) [@jpfeifer](https://forum.dynare.org/u/jpfeifer)\
**Post date:** [11 April 2014 13:26 UTC](https://forum.dynare.org/t/simulated-moments/670/11 "2014-04-11T13:26:42Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![quiko](https://forum.dynare.org/letter_avatar_proxy/v4/letter/q/49beb7/32.png) [@quiko](https://forum.dynare.org/u/quiko)\
**Post date:** [11 April 2014 14:11 UTC](https://forum.dynare.org/t/simulated-moments/670/12 "2014-04-11T14:11:17Z")

</div>

Thanks a lot, it’s perfect.

---

<div class="post-metadata">

**Author:** ![quiko](https://forum.dynare.org/letter_avatar_proxy/v4/letter/q/49beb7/32.png) [@quiko](https://forum.dynare.org/u/quiko)\
**Post date:** [2 July 2014 21:48 UTC](https://forum.dynare.org/t/simulated-moments/670/13 "2014-07-02T21:48:02Z")

</div>

[quote=“jpfeifer”]It depends what exactly you need. Variances are stored in oo\_.variance, so

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

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.

---

<div class="post-metadata">

**Author:** ![jpfeifer](https://forum.dynare.org/user_avatar/forum.dynare.org/jpfeifer/32/5044_2.png) [@jpfeifer](https://forum.dynare.org/u/jpfeifer)\
**Post date:** [3 July 2014 09:58 UTC](https://forum.dynare.org/t/simulated-moments/670/14 "2014-07-03T09:58:06Z")

</div>

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

```auto

```

---

<div class="post-metadata">

**Author:** ![quiko](https://forum.dynare.org/letter_avatar_proxy/v4/letter/q/49beb7/32.png) [@quiko](https://forum.dynare.org/u/quiko)\
**Post date:** [3 July 2014 11:19 UTC](https://forum.dynare.org/t/simulated-moments/670/15 "2014-07-03T11:19:40Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![quiko](https://forum.dynare.org/letter_avatar_proxy/v4/letter/q/49beb7/32.png) [@quiko](https://forum.dynare.org/u/quiko)\
**Post date:** [3 July 2014 12:05 UTC](https://forum.dynare.org/t/simulated-moments/670/16 "2014-07-03T12:05:30Z")

</div>

[quote=“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?: covars(:,s)=(oo\_.var(1,:));

Thank you for your replies.[/quote]

---

<div class="post-metadata">

**Author:** ![jpfeifer](https://forum.dynare.org/user_avatar/forum.dynare.org/jpfeifer/32/5044_2.png) [@jpfeifer](https://forum.dynare.org/u/jpfeifer)\
**Post date:** [4 July 2014 08:37 UTC](https://forum.dynare.org/t/simulated-moments/670/17 "2014-07-04T08:37:06Z")

</div>

Then you are looking for something like

`temp=triu(oo_.var,1)
covars(:,s)=temp(temp~=0);`  
Note that the vectorization is in row major order, i.e. you get the covariance of 1 with 2 first, then 1 with 3, 2 with 3, 1 with 4, 2 with 4 etc.

---

<div class="post-metadata">

**Author:** ![hlt215](https://forum.dynare.org/letter_avatar_proxy/v4/letter/h/3bc359/32.png) [@hlt215](https://forum.dynare.org/u/hlt215)\
**Post date:** [20 December 2014 22:34 UTC](https://forum.dynare.org/t/simulated-moments/670/18 "2014-12-20T22:34:22Z")

</div>

Hi again,

A further follow-up question about this. Instead of one loop over a parameter, I now have an additional loop over another parameter. Thus, two different loops in the same simulation. So I have:

gammas = 0.3:0.1:1;  
for s = 1:length(gammas);  
gamma = gammas(s);

alphas = 0.3:01:1;  
for z = 1:length(alphas);  
alpha = alphas(z);

stoch\_simul(hp\_filter = 1600, order = 1, nograph);

end;  
end;

Previously, I extracted the variance by doing:  
moments(:,s)=(diag(oo\_.var));

How should I change the above command to account also for the other loop? I tried to do like this:  
moments(:,t)=(diag(oo\_.var)); where t = s_z. But I would like to put my results into a matrix of dimension s_z in which each element is the variance evaluated at a the specific pair of (s,z).

---

<div class="post-metadata">

**Author:** ![jpfeifer](https://forum.dynare.org/user_avatar/forum.dynare.org/jpfeifer/32/5044_2.png) [@jpfeifer](https://forum.dynare.org/u/jpfeifer)\
**Post date:** [28 December 2014 10:52 UTC](https://forum.dynare.org/t/simulated-moments/670/19 "2014-12-28T10:52:23Z")

</div>

The variances are a vector (1 dimension). Now you want to store this vector along two additional dimensions. The result will be three dimensional. How can you store it in a two dimensional matrix? This would only be possible if the variance is a scalar.

---

<div class="post-metadata">

**Author:** ![hlt215](https://forum.dynare.org/letter_avatar_proxy/v4/letter/h/3bc359/32.png) [@hlt215](https://forum.dynare.org/u/hlt215)\
**Post date:** [28 December 2014 13:15 UTC](https://forum.dynare.org/t/simulated-moments/670/20 "2014-12-28T13:15:47Z")

</div>

Thanks. What about if I only look at one variable? For example, the command moments(:,t) gives me all variables in the rows. If I look only at one variable, should I be able to get the variance stored in a matrix of dimension (s,z)? That is, each element of the matrix to be the variance evaluated at a specific pair (s,z)?  
I exactly need a three dimension matrix because I would like to make a 3D plot. (I want to see how the variance changes when both parameters change at the same time).

[Next page](https://forum.dynare.org/t/simulated-moments/670.md?page=2)