Attempting ARMA(1,1) noise regression model

[Wayne Folta]

I'm calculating the trend of a climate temperature time series, which will have serial correlation, so I've tried to create a regression which has ARMA(1,1) noise, which I've also heard is basically equivalent to AR(1) noise + white noise. My model (in which I think I'm creating AR(1) + white noise) is:

data
{
int<lower=1> N ;
real x[N] ;
real y[N] ;
}

parameters
{
real alpha ;
real beta ;
real kappa ;

real<lower=0> sigma0[N+1] ;
real<lower=0> sigma1 ;
}

model
{
real sigma ;

alpha ~ cauchy (0, 5) ;
beta ~ cauchy (0, 5) ;

kappa ~ gamma (1.2, 1) ;
sigma1 ~ gamma (3, 1) ;

for (n in 1:(N+1)) { sigma0[n] ~ gamma (3, 1) ; }

sigma <- sigma1 ;

for (n in 1:N)
{
sigma <- fabs (sigma0[n] + kappa * sigma) ;

y[n] ~ normal (alpha + beta * x[n], sigma) ;
}
}

I've tried several variations on this, but the bottom line is I'm not sure how to handle the scale parameter, sigma, properly. There's an example in the manual of an AR(1) series, but not of AR(1) noise, and I'm stumped. Also, I don't think the fabs is a good idea, but in various versions of the code I've ended up with errors where sigma was negative. Any help appreciated!

(I'm comparing my results to using R's `arima` with an xreg, and they definitely differ, though if I throw ablines onto a graph of the data, each looks "reasonable".)

Thanks!

Wayne

[Bob Carpenter]

:-) I just added a MA model example to Stan yesterday.
There's a hard-coded MA(2) example and a general
MA(Q) example included:

https://github.com/stan-dev/stan/tree/master/src/models/misc/moving-avg

For a start, your model can be full vectorized, which will run much
faster. For instance, this:

> for (n in 1:(N+1)) { sigma0[n] ~ gamma (3, 1) ; }

can be coded as

sigma0 ~ gamma(3,1);

And if you declare x as

vector[N] x;

then you can replace this:

> for (n in 1:N)

> y[n] ~ normal (alpha + beta * x[n], sigma) ;

with this more efficient equivlaent

y ~ normal(alpha + beta * x, sigma);

I'm afraid I'm not sure what your intent is with this:

> sigma <- fabs (sigma0[n] + kappa * sigma) ;
>
> y[n] ~ normal (alpha + beta * x[n], sigma) ;

You've constrained sigma0 and sigma to be positive, so the
result is only guaranteed to be positive if you constrain kappa
to be positve, too:

real<lower=0> kappa;

Absolute values aren't ideal because of the non-monotonic
effect on sigma; it's more usual to use a log-linear model
for positive values, setting

sigma <- exp(...);

- Bob

[Wayne]

Thanks, I've implemented what you've suggested. At least the parts I understand. It's now usually close to what I get with R's arima (with xreg) now, but then there are times where it's obviously different.

And as I've thought about it, I think the problem is that the normal distribution is symmetric, so whatever I do with the scale will increase the spread but won't push it in one direction or another. Which an AR(1) noise would do, I think: tend to push things in one direction for a while, then another (positive AR coefficient), or back and forth (negative AR coefficient). So whether I constrain kappa or use fabs or whatever, the AR(1) noise is too free to move "the wrong way". At least that's how I think it works.

So my guess is that the model works well when the AR coefficient it positive and not so well when it's negative. I've thought about using a beta distribution, scaled appropriately to extend outside of [0,1], but it's pretty complicated to imagine how to make everything work appropriately.

[Wayne]

Bob,

I believe I have a model working with an AR(1) noise:

data

    {

    int<lower=1> N ;

    real x[N] ;

    real y[N] ;

    }

parameters

    {

    real alpha ;

    real beta ;

    real kappa ;

    real lambda ;

    real<lower=0> sigma ;

    }

model

    {

    real error ;

    alpha ~ normal (0, 30) ;

    beta  ~ normal (0, 2) ;

    kappa ~ normal (0, 2) ;

    lambda ~ normal (0, 2) ;

    sigma ~ gamma (3, 1) ;

    error <- 0 ;

    for (n in 1:N)

{

y[n] ~ normal (alpha + beta * x[n] + error, sigma) ;

error <- kappa * (y[n] - (alpha + beta * x[n])) ;

}

    }

[Bob Carpenter]

Cool. Can I borrow this for the manual, with attribution
to you, of course? I ran out of steam before writing down
a full ARMA model. Even better if you have an R program ot
to simulate data or a public data set with a known answer
I can include in the distro.

I really like the way you defined the error:

>
> error <- 0 ;
>
> for (n in 1:N) {
>
> y[n] ~ normal (alpha + beta * x[n] + error, sigma) ;
>
> error <- kappa * (y[n] - (alpha + beta * x[n])) ;

I'm finding Stan user programs' use of curly braces in blocks
and space in expressions like a new art form!

- Bob

> --
> You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Wayne]

Bob,

I posted something a bit less-than-helpful a couple of minutes ago and attempted to delete it, but wasn't able to. Sorry for the spam!

I created a test case:

sa <- arima.sim (list (ar=c(0.7), order=c(1, 0, 0)), 1000) + seq (0.5, 500, 0.5)
plot (sa)
stan.list1 <- list (N=length (sa), x=1:1000, y=sa)


stan.model1 <- stan (model_code=stan.code1, model_name="sa", data=stan.list1, iter=20000, chain=4)


print (stan.model1, digits_summary=8)

Which has a slope of 0.5 and an intercept of 0, with AR(1) noise and an AR coefficient of 0.7. With STAN, I got reasonable results:

              mean    se_mean         sd 

alpha   -0.1293778 0.00181275 0.21566235 

beta     0.5004723 0.00000317 0.00037522 

kappa    0.6980644 0.00021365 0.02354209 

The Rhat's were all less than 1.0002 which is a good sign as well.

In case the system did manage to delete my previous posting (which was basically the same, except with a regression intercept and slope of zero which is not very visible), you have my permission to use my example code however you want. I'll note that I'm still learning all of this stuff, so passing it by a statistical expert would be a good idea.

[Bob Carpenter]

Nice. Thanks!

- Bob

On 3/11/13 10:51 PM, Wayne wrote:
> I just tried generating my own test case:

>
> |
> sa <- arima.sim (list (ar=c(0.7), order=c(1, 0, 0)), 1000)
> |
>

> and get the official answer:
>
> |
> arima (sa, c(1, 0, 0), xreg=1:1000)
> |
>
> which yielded:
>
> ar1 intercept 1:1000
> 0.7164 0.0346 -1e-04
> s.e. 0.0220 0.2126 4e-04
>
>
> and what I got from STAN was:
>
> mean se_mean sd
> alpha 0.01980100 0.00171910 0.21660766
> beta -0.00008250 0.00000293 0.00037689
> kappa 0.72007884 0.00018809 0.02241440
>
> Which is pretty close. I haven't tested it exhaustively, and I'd recommend that before putting it into the manual, but I
> believe it's right, and you certainly have my permission. This is a great tool.

>
> On Monday, March 11, 2013 8:56:03 PM UTC-4, Bob Carpenter wrote:
>

> Cool. Can I borrow this for the manual, with attribution
> to you, of course? I ran out of steam before writing down
> a full ARMA model. Even better if you have an R program ot
> to simulate data or a public data set with a known answer
> I can include in the distro.
>