[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

From |
Marcello Pagano <[email protected]> |

To |
[email protected] |

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" <[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

hassomewhat 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 dois: ********************* 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)andthe rest as epsilons */Using the "burn-in" method Sasha just described, just use the first `k'

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?

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/

* * 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/

**References**:**Re: st: time series simulation study***From:*[email protected] (Jeff Pitblado, Stata Corp.)

**Re: st: time series simulation study***From:*"Sascha O. Becker" <[email protected]>

- Prev by Date:
**RE: st: RE: bug with histogram?** - Next by Date:
**RE: st: RE: Weibull simulation-parametrization** - Previous by thread:
**Re: st: time series simulation study** - Next by thread:
**st: truncation and endogenous stratification in TCM** - Index(es):

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