Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Re: st: time series simulation study


From   "Sascha O. Becker" <[email protected]>
To   [email protected]
Subject   Re: st: time series simulation study
Date   Wed, 7 May 2003 10:43:56 +0200 (MEST)

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/



© Copyright 1996–2024 StataCorp LLC   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index