# Re: st: time series simulation study

 From Marcello Pagano To statalist@hsphsun2.harvard.edu Subject Re: st: time series simulation study Date Wed, 07 May 2003 08:31:11 -0400

I don't know if the question makes sense, but the answer makes even less sense. y(t) as defined by Sascha is not an autoregression. It is not stationary (requirement for "stabilisation") for any set of b (except for b0=0). The solution does not satisfy the model definition, as given by Sascha. The solution attempts to generate (I say attempts because this dependence on the dying out of the starting conditions is a pretty sloppy way to proceed, especially if k is large and any of the roots of the generating polynomial for the bs are anywhere near the unit circle,
when a much better solution is to generate the first k values according to a multivariate distribution with the proper covariance matrix, and then solve the difference equation.)
y(t) = mean(t) + signal(t)
where the signal is autoregressive.
This is not what satisfies Sascha's equation. Maybe this is what Sascha wants?

m.p.

Sascha O. Becker wrote:

```Thanks a lot , Jeff!!

```
```"Sascha O. Becker" <so.b@gmx.net> asks how to simulate an AR(k) process
when
k>1:

```
Dear statalisters,

in a simulation study (with something like 1000 replications), I want to
generate simple AR(k) time series
y(t) = c + b0 * t + b1 * y(t-1) + b2* y(t-2) +...+ bk* y(t-k)+ eps.

where I want to pick values t=450, ..., 500, i.e. after the time series

```has

```
somewhat stabilised. I.e. I want to generate 1000 time series like this
(using simul "around" the procedure below) .
For the case where b2=...=bk=0, i.e. an AR(1), I can do it pretty easily
using the feature that Stata's "replace" command replaces recursively.

```What I do

```
```is:

*********************
set seed 1;
set obs 10;
generate mynormal=invnorm(uniform());
g time=_n;
tsset time;
g y=.;
replace y=mynormal in 1; /* use the first entry in mynormal as a y(0)

```
```and

```
the rest as epsilons */
replace y=1+b0*time+b1*l.y+mynormal if time>1; /* NOTE: replace works
recursively !! */
list;
*********************

which yields

list;

mynormal time y
1. .4349075 1 .4349075
2. 1.460773 2 13.0046
3. -.0389092 3 29.03223
4. -.0759261 4 49.95335
5. 1.921223 5 77.76375
6. .1418411 6 108.6329
7. .8990735 7 145.0734
8. .9799716 8 186.3653
9. .2037313 9 231.6133
10. -.3139055 10 281.0451

I do not see immediately how to generate higher order AR(k) (k>1). Any suggestions?

Using the "burn-in" method Sasha just described, just use the first `k'
normal
errors as the starting values (this was the procedure for `k'==1). I made
some modifications to the above snippet of code to do just that. Notice
that
I generate the zero mean AR(3) process first, then add in the trend.

***** BEGIN
* sim-ar.do

version 7

clear
set mem 10m
set seed 1
local keep 200
local burn 50000
local obs = `keep' + `burn'
set obs `obs'
* error term
generate double mynormal = invnorm(uniform())
* time variable
generate double time = _n - `burn'
tsset time
* AR(3) coefficient vector, with iid normal noise
matrix b_ar = .2, -1, .1, 1
matrix colnames b_ar = l.y l2.y l3.y mynormal
local k = colsof(b_ar)
local kp1 = `k'+1
* generate time-series response
generate double y = mynormal in 1/`k'
matrix score y = b_ar in `kp1'/l, replace
* keep "burned-in" observations
keep in -200/l
* add the other regression structure last
matrix b0 = 3, -1
matrix colnames b0 = time _cons
matrix score double xb = b0
replace y = y + xb
* check coefficients using arima
arima y time , ar(1/3)

exit
***** END

Here are the results from running the above do-file (I ran it using both
Stata
7 and Stata 8).

***** BEGIN
. * sim-ar.do
. . version 7

. . clear

. set mem 10m
(10240k)

. set seed 1

. local keep 200

. local burn 50000

. local obs = `keep' + `burn'

. set obs `obs'
obs was 0, now 50200

. * error term
. generate double mynormal = invnorm(uniform())

. * time variable
. generate double time = _n - `burn'

. tsset time
time variable: time, -49999 to 200

. * AR(3) coefficient vector, with iid normal noise
. matrix b_ar = .2, -1, .1, 1

. matrix colnames b_ar = l.y l2.y l3.y mynormal

. local k = colsof(b_ar)

. local kp1 = `k'+1

. * generate time-series response
. generate double y = mynormal in 1/`k'
(50196 missing values generated)

. matrix score y = b_ar in `kp1'/l, replace

. * keep "burned-in" observations
. keep in -200/l
(50000 observations deleted)

. * add the other regression structure last
. matrix b0 = 3, -1

. matrix colnames b0 = time _cons

. matrix score double xb = b0

. replace y = y + xb
(200 real changes made)

. * check coefficients using arima
. arima y time , ar(1/3)

(setting optimization to BHHH)
Iteration 0: log likelihood = -305.2692 Iteration 1: log likelihood = -304.35313 Iteration 2: log likelihood = -304.2739 Iteration 3: log likelihood = -304.26906 Iteration 4: log likelihood = -304.26866 (switching optimization to BFGS)
Iteration 5: log likelihood = -304.26855 Iteration 6: log likelihood = -304.26851 Iteration 7: log likelihood = -304.26851
ARIMA regression

Sample: 1 to 200 Number of obs = 200
Wald chi2(4) = 1.36e+07
Log likelihood = -304.2685 Prob > chi2 = 0.0000

```------------------------------------------------------------------------------

```
```            |                 OPG
y            |      Coef.   Std. Err.      z    P>|z|     [95% Conf.
Interval]

```
```-------------+----------------------------------------------------------------

```
y |
time | 2.999811 .0008205 3656.29 0.000 2.998203 3.001419
_cons | -.9472167 .0992725 -9.54 0.000 -1.141787 -.7526462

```-------------+----------------------------------------------------------------

```
ARMA |
ar |
L1 | .1905617 .0727386 2.62 0.009 .0479967 .3331267
L2 | -1.000933 .0092656 -108.03 0.000 -1.019093 -.9827729
L3 | .0942806 .0728187 1.29 0.195 -.0484414 .2370025

```-------------+----------------------------------------------------------------

```
/sigma | 1.085205 .0583985 18.58 0.000 .9707459 1.199664

```------------------------------------------------------------------------------

```
. . exit
***** END

--Jeff
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/

```
*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/
```