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/
>
--
+++ GMX - Mail, Messaging & more http://www.gmx.net +++
Bitte l�cheln! Fotogalerie online mit GMX ohne eigene Homepage!
*
* 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/
