Specify 'autocor' as a term inside the formula

[paul-buerkner]

Currently, the autocor argument is the only place outside of formula where R formulas are used to express parts of the model. This has historical reasons, but I am no longer sure if keeping autocor separate is actually beneficial. As an alternative, it might make sense to define a new special function, ac say, that behaves like other special terms such as mo or me but handles autocorrelation structures. For some most relevant structures, such a CAR and some ARMA versions, this would actually reflect much better what actually happens inside, which is simply adding a term to the linear predictor. for an ARMA(1, 1) model, this could look for instance as follows:

y ~ ac(time, type = "arma", p = 1, q = 1) 

and we could also add convient wrappers around ac such as

y ~ arma(time, p = 1, q = 1) 

Not sure yet about all the internal implications, but it may be an idea worth thinking about further.

[ASKurz]

This has a nice intuitive feel, to me

[SteveBronder]

I like this a lot! For distributional parameters this becomes fun if there can be something like

sigma ~ garch(time | group, 1, 1)

[paul-buerkner]

Incidentally, I started working on this feature yesterday and it is good to have someone to discuss this with.

I agree it is intuitive to specify residual related autocors in the formula of sigma, but in terms of implementation, this is not entirely what we want. The reason is that we do not really 'predict' sigma, but either add the residual predictors to the mean arma(..., cov = FALSE) or model the residual correlation structure arma(..., cov = TRUE) neither of which actually predicts sigma which, as a parameter, still remains constant across observations.

[SteveBronder]

np It's a fun problem!

It feels like I sort of cheated with the above since garch is made pretty for forecasting on sigma and would probably gen something like the below for cov=FALSE

  for (n in 1:N) {
   // standard arma stuff
    mu[n] += Err[n, 1:Kma] * ma;
    err[n] = Y[n] - mu[n];
    for (i in 1:J_lag[n]) {
      Err[n + 1, i] = err[n + 1 - i];
    }
    mu[n] += Err[n, 1:Kar] * ar;
    // garch stuff
    // pretend we have a pow that does elementwise
    sigma[n] += pow(Err[n, 1:GKar], 2) * g_ar;  
    for (i in 1:G_lag[n]) {
      sigma_lag[n + 1, i] = sigma[n + 1 - i];
    }
    sigma[n] += pow(sigma_lag[n, 1:GKrr], 2) * g_rr;
  }

What other forms instead of ar() were you thinking about? I could also see a separate block of formulas for modeling residuals

or model the residual correlation structure arma(..., cov = TRUE)

I have to look more into the residual correlation structure code before I can say anything meaningful about how to transfer that from user space to gen code.

Maybe @asael697 would have some opinions on this since he wrote the varstan package

[paul-buerkner]

Yes, GARCH indeed looks like a perfect application for predictions of sigma. I wasn't sure what exactly GARCH does but now it makes a lot of sense. Thank you.

[SteveBronder]

Could it be broken down by parameter type? Like continuous unbounded parameters whose location and scale are univariate derivations from the mean and variance have a certain set of generic lag functions like ar, ma, (diff?), garch and it's wackier variants.

[asael697]

Hi in my package varstan, I contemplate the seasonal arima and garch classes where you can do:

• Dynamic regressions with arima errors.
  -garch
• arma-garch
• Dynamic harmonic regressions
• Seasonal arima
  The official version 1 will be online soon(this days).

And I am currently working on the wackier variants for garch models

[asael697]

Saddly I dont contemplate multi-level time series models, or stochastic volatility models. But maybe in my version 2 I will update some more nice stuff

[SteveBronder]

Nice and no worries!

[asael697]

np It's a fun problem!

It feels like I sort of cheated with the above since garch is made pretty for forecasting on sigma and would probably gen something like the below for cov=FALSE

  for (n in 1:N) {
   // standard arma stuff
    mu[n] += Err[n, 1:Kma] * ma;
    err[n] = Y[n] - mu[n];
    for (i in 1:J_lag[n]) {
      Err[n + 1, i] = err[n + 1 - i];
    }
    mu[n] += Err[n, 1:Kar] * ar;
    // garch stuff
    // pretend we have a pow that does elementwise
    sigma[n] += pow(Err[n, 1:GKar], 2) * g_ar;  
    for (i in 1:G_lag[n]) {
      sigma_lag[n + 1, i] = sigma[n + 1 - i];
    }
    sigma[n] += pow(sigma_lag[n, 1:GKrr], 2) * g_rr;
  }

What other forms instead of ar() were you thinking about? I could also see a separate block of formulas for modeling residuals

or model the residual correlation structure arma(..., cov = TRUE)

I have to look more into the residual correlation structure code before I can say anything meaningful about how to transfer that from user space to gen code.

Maybe @asael697 would have some opinions on this since he wrote the varstan package

Returning to the garch issues I also think add a garch structure is a good idea, even so garch is for predictions on sigma, the principal goal of garch's and SVM is to model the heteroscedasticity of the residuals, property that happens when the square residuals have "large" autocorrelation.

[asael697]

Could it be broken down by parameter type? Like continuous unbounded parameters whose location and scale are univariate derivations from the mean and variance have a certain set of generic lag functions like ar, ma, (diff?), garch and it's wackier variants.

One of the big problems that I have with differencing the time series, is that:

• It transforms your data,
• the predictions with differenced data are really bad, and gets worse if you increment the number of diff (and that is why nobody likes arima models).

At the begining I try to solve it inside stan, differencing the data in stan and apply the arma model. But it give me alot of divergencies. So now I just give the transformed data to stan. (Also you have to implement diff and diffinv functions inside stan, and I still think my version were inefficients ).

Other big problem with diff and dynamic regression is that you have to apply diff to y (dependent) and to x (the regressor variables)

[SteveBronder]

Returning to the garch issues I also think add a garch structure is a good idea, even so garch is for predictions on sigma, the principal goal of garch's and SVM is to model the heteroscedasticity of the residuals, property that happens when the square residuals have "large" autocorrelation.

Would be v nice to add! I think it could fit in the current model, though it would be fun to have arma on the mean and a garch on sigma. idk how that would work currently tho.

There are AR schemes for poisson, gamma, and a few other distributions. Though I think those would be out of scope for an initial scheme. Could be a good GSOC thing

[asael697]

There are also for garch, the intgarch (integer-garch) process.
@SteveBronder One of my fears of the wackier variants are the e-garch and aparch models, because of the aboslute value, and is an absolute value on a errorr(0,1) close to zero, so that could bring alot of divergencies on the Hamiltonian montecarlo (at the moment of evaluating the jacobian). Still I haven´t try it yet.

For the e-garch is not big deal, cause instead of an e-garch(exponential garch) you could use an SVM (stochastic volatility model) that is similar and have a better performance.

[SteveBronder]

There are also for garch, the intgarch (integer-garch) process.

Nice!

@SteveBronder One of my fears of the wackier variants are the e-garch and aparch models, because of the aboslute value, and is an absolute value on a errorr(0,1) close to zero, so that could bring alot of divergencies on the Hamiltonian montecarlo (at the moment of evaluating the jacobian). Still I haven´t try it yet.

Yeah I could see that being bad news bears. I'm also not sure exactly what the egarch would look like. In the link below is z_t supposed to be the conditional error? Or would we literally do

matrix[N, Sims] z = std_normal_rng()

https://vlab.stern.nyu.edu/docs/volatility/EGARCH

[asael697]

I think that should work, cause z_t is a white noise so actually it doesn't have to be conditioned to the previous error lags

[paul-buerkner]

I have just merged PR #832 that implements this feature and so I will close this issue. @SteveBronder would you mind opening a new issue for the GARCH model and summarize your state of discussion about these models there? I kind of lost you both in the middle of it.

[SteveBronder]

Awesome! Yeah sure I'll open up a new issue for the garch stuff