Thanks a lot , Jeff!!
"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
(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
[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/
