# Re: st: time series simulation study

 From [email protected] (Jeff Pitblado, Stata Corp.) To [email protected] Subject Re: st: time series simulation study Date Tue, 06 May 2003 10:37:24 -0500

```"Sascha O. Becker" <[email protected]> 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

. * 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
[email protected]
*
*   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/
```