HMM and ARCH(1) model

[Wad Wan]

Hello guys,

I want to do ARCH(1) model for 2-state Markov Regime Switching Model for some financial time series.

In Stan Modeling Language User's Guide and Reference Manual (v.2.6.2) on page 69
I have found how to implement ARCH(1) model. This code works well, and I get correct approximations.

On page 81 I find the HMM "Semisupervised Estimation", but I can't figure out how to combine the ARCH(1)
on page 61 and Semisupervised Estimation on p. 81.

Does anyone have an idea how to do it?..

All the best, Serge

[Wad Wan]

Yo-yo, Stan, check my rhyme!..

So, guys to be more specific, I mean the following. Suppose that we have a financial or economic time series, where we have a prior subjective belief that there are two different latent states (i.e., bull and bear markets). Now, we want to use a hidden Markov models (i.e., Markov Regime State Switching) that is somewhat similar to SWARCH model by Hamilton and Susmel. But to keep it simple and stupid, we begin to look at ARCH(1) model and Markov Regime State Switching Model. Now, what I have a great difficulty to figure out how to use Stan to combine these two models together. When I read Stan's user manual, they write very well the explanation for how to solve ARCH(1) model, and the explanation for how to solve HMM, but there is that missing link that how do you combine these two models in Stan, such that, we run HMM recursively through the time series, and use ARCH(1) model to optimize parameters in the two latent states?...

So, to keep it simple and stupid we begin with the first order ARCH(1) stochastic process estimation of paramters in the model. The definition of ARCH(1) looks something like that:

\begin{eqnarray}
     a_t = \sigma_t \epsilon_t, \text{ } \sigma_t = \sqrt[]{\alpha_0 + \alpha_1 a^2_{t-1}},  \text{ }  \sigma_t^2 = \alpha_0 + \alpha_1 a^2_{t-1},
 \end{eqnarray}

where $\alpha_0 > 0$, $\alpha_1 \geq 0$, $\epsilon_t$ is independent of the past, and $\sigma_t$ is known function value of the prior information
at time $t-1$. 

We begin our implementation in R and Stan, which looks like that:

######################################################################
## 1st we simulate ARCH(1) stochastic process, with known user-defined-true-parameters
######################################################################

  library(rstan)   # we load in the Stan package
  library(tseries)

  alpha0 <- 0.06321   # user defined subjective constant
  alpha1 <- 0.03548   # user defined subjective constant
  T <- 10000               # number of replications

  r_t <- rep(0, T)

  sigma2_t <- 0.1        # initial guess; user defined parameter
  r_t[1] <- sqrt(sigma2_t)*rnorm(1,0,1)

  for (i in 2:T) {
     eps <- rnorm(1,0,1)   
     sigma2_t <- alpha0 + alpha1*r_t[i-1]^2
     r_t[i] <- sqrt(sigma2_t)*eps
  }

  plot.ts(r_t)    # vizualization of the ARCH(1) process
  hist(r_t)


  ## 2nd step: Estimation of ARCH(1) model's latent paramters, by using Stan software:
  ######################################################################
  ## Stan modelling in R
  ######################################################################

  modelString = "

     data {
         int<lower=0> T; // number of time points
         real r_t[T]; // return at time t
     }

     parameters {
         real mu;
         real<lower=0> alpha0; // noise intercept
         real<lower=0,upper=1> alpha1; // noise slope
      }

    model {

        for (t in 2:T)
            r_t[t] ~ normal(mu, sqrt(alpha0 + alpha1 * pow(r_t[t-1] - mu,2)));
    }
 "

  sm <- stan_model(model_code = modelString) # Stan model
  fit <- sampling(sm, data=list(r_t=r_t, T=T), iter=1000, chains=4, warmup=250, thin=5); fit
  print(fit, digits_summary = 5)
  plot(fit)

  fit_sim <- extract(fit)
  hist(fit_sim$alpha0)
  hist(fit_sim$alpha1)

  ######################################################################
  ## Done!..
  ######################################################################

Sub-conclusion: Stan estimates the latent variables of ARCH(1) model, i.e., 'alpha0' and 'alpha1'
correctly, although there might be slight decimal underestimation of the 'alpha1' parameters.

So, the second step is to use hidden Markov models (i.e. Markov Regime State Switching) to find
the hidden states. As already mentioned above, there is a nice explanation example in Stan user's manual,
but it is not exhastive with ARCH(1) modelling. Now, I have also tried to look at the following link
https://github.com/jrnold/ssmodels-in-stan/blob/master/stan/hmm.stan by Jeffrey Arnold, but could't really
figure out how to combine ARCH(1) model with HMM latent state estimation. His code looks like that

/*
   Hidden Markov Model Example

   The following model

   Let $p(\theta_1)$, and for $t = 2, \dots ,T$,
   $$
   \begin{aligned}[t]
   p(\theta_t | \theta_{1:t-1}, \psi) = p(\theta_t | \theta_{t-1}, \psi)  \\
   p(y_t | y_{1:t-1}, \theta_{1:t}) = p(y_t | \theta_t)
   \end{aligned}
   $$

   Stan does not allow sampling discrete parameters.
   Thus, the sampling procedes in two steps,

   1. In the ``model`` block, the observed likelihood is calculated by marginalizing out the discrete latent states.
   2. In the ``generated quantities`` block, the discrete states are sampled conditional on the estimated parameters.  Additionally, statistics of the discrete states are calculated: the conditional probabilities of the states decoding $\Pr(\theta_{t} | y, \psi)$, the global decoding (Viterbi algorithm) $\argmax p(\theta_{1:T} | y_{1:T}, \psi)$.
  
 */
data {
  int n; // time
  int m; // number of states
  int y[n];
}
parameters {
  // transition probabilities
  simplex[m] Gamma[m]; // Gamma[from, to]
  // probability parameters
  positive_ordered[m] mu;
  real<lower = 0.0> mu_a;
  real<lower = 0.0> mu_b; 
}
transformed parameters {
  vector<upper = 0>[m] logp[n];
  real llik;
 
  simplex[m] delta;
  matrix[m, m] Gamma_mat;
  for (i in 1:m) {
    for (j in 1:m) {
      Gamma_mat[i, j] <- Gamma[i, j];
    }
  }
 
  for (i in 1:n) {
    for (j in 1:m) {
      logp[i, j] <- poisson_log(y[i], mu[j]);
    }
  }
  // calculate the stationary initial distribution
  delta <- (rep_row_vector(1.0, m) /
          (diag_matrix(rep_vector(1.0, m))
        - Gamma_mat + rep_matrix(1.0, m, m))) ';
 
  // HMM Log likelihood
  // Lystig and Hughes (2002) "Exact Computation ..."
  {
    matrix[n, m] log_lambda;
    vector[n] log_Lambda;
    for (j in 1:m) {
      log_lambda[1, j] <- log(delta[j]) + logp[1, j];
    }
    log_Lambda[1] <- log_sum_exp(log_lambda[1]);
    for (t in 2:n) {
      for (j in 1:m) {
    vector[m] tmp;
    for (i in 1:m) {
      tmp[i] <- (log_lambda[t - 1, i]
             + log(Gamma[i, j]));
    }
    log_lambda[t, j] <- (log_sum_exp(tmp)
                 + logp[t, j]
                 - log_Lambda[t - 1]);
      }
      log_Lambda[t] <- log_sum_exp(log_lambda[t]);
    }
    llik <- sum(log_Lambda);
  }
}
model {
  increment_log_prob(llik);
  // regularize the states.
  mu ~ gamma(mu_a, mu_b);
}
generated quantities {
  // log forward probabilities
  vector[m] log_alpha[n];
  // log log backward probabilities
  vector[m] log_beta[n];
  // state probabilities
  vector[m] stateprob[n];
  // global decoding, calcualated with Viterbi algorithm
  int viterbi[n];
  // sample of latent states
  int states[n];

  // Forward probabilities
  // this could be made more efficient by reusing forward probabilities
  // calculated in the likelihood
  {
    real u;
    vector[m] phi;
    vector[m] v;
    real lscale;
    phi <- delta;
    lscale <- 0.0;
    for (i in 1:n) {
      if (i == 1) {
          v <- exp(log(delta) + logp[1]);
      } else {
      v <- exp(log(phi ' * Gamma_mat) + logp[i] ') ';
      }
      u <- sum(v);
      lscale <- lscale + log(u);
      phi <- v / u;
      log_alpha[i] <- log(phi) + lscale;
    }
  }
  // backward probilities
  {
    int t;
    real u;
    vector[m] phi;
    vector[m] v;
    real lscale;

    log_beta[n] <- rep_vector(0.0, m);
    phi <- rep_vector(1.0 / m, m);
    lscale <- log(m);
    for (i in 1:(n - 1)) {
      t <- n - i;
      v <- Gamma_mat * exp(logp[t + 1] + log(phi));
      log_beta[t] <- log(v) + lscale;
      u <- sum(v);
      phi <- v / u;
      lscale <- lscale + log(u);
    }
  }
  // state conditional probabilities
  // uses forward and backward probabilities
  for (i in 1:n) {
    stateprob[i] <- exp(log_alpha[i] + log_beta[i] - llik);
  }
  // Sample values from the states
  states[n] <- categorical_rng(exp(log_alpha[n]) / sum(exp(log_alpha[n])));
  for (i in 1:(n - 1)) {
    int t;
    vector[m] theta;
    t <- n - i;
    for (j in 1:m) {
      theta[j] <- exp(log_alpha[t, j]) * Gamma[j, states[t + 1]];
    }
    theta <- theta / sum(theta);
    states[t] <- categorical_rng(theta);
  }
  // Global Decoding (Viterbi algorithm)
  //
  {
    vector[m] log_xi[n];
    // forwards pass
    log_xi[1] <- log(delta) + logp[1];
    for (t in 2:n) {
      for (j in 1:m) {
    real max_xi_gamma;
    real tmp;
    max_xi_gamma <- negative_infinity();
    for (i in 1:m) {
      tmp <- log_xi[t - 1, i] + log(Gamma[i, j]);
      if (tmp > max_xi_gamma) {
        max_xi_gamma <- tmp;
      }
    }
    log_xi[t, j] <- max_xi_gamma + logp[t, j];
      }
    }
    // backwards pass
    {
      real tmp;
      tmp <- negative_infinity();
      for (i in 1:m) {
        if (log_xi[n, i] > tmp) {
      tmp <- log_xi[n, i];
          viterbi[n] <- i;
        }
      }
      for (s in 1:(n - 1)) {
        int t;
        real tmp1;
        t <- n - s;
        tmp1 <- negative_infinity();
        for (i in 1:m) {
          real tmp2;    
          tmp2 <- log_xi[t, i] + log(Gamma[i, viterbi[t + 1]]);
          if (tmp2 > tmp1) {
        tmp1 <- tmp2;
            viterbi[t] <- i;
          }
        }
      }
    }
  }
}

Now, here come the question. How do you take the code for ARCH(1) model and combine it with the HMM code (e.g., Jeffrey Arnold's code). I mean more like,
when HMM recursively changes between states, we need to run ARCH(1) for every state?... Right?... Or, correct me if I'm wrong, did I misunderstood it all?..

All the best, Serge Boris Nesterenko

[Bob Carpenter]

Sorry this question has languished for so long. I still
don't quite understand the intent, so I haven't answered.

> we run HMM recursively through the time series, and use
> ARCH(1) model to optimize parameters in the two latent states?...

I wrote the examples for the manual, but am hardly a time-series
expert. Usually if you can write the model's likelihood down, it's
down to some programming. Which can be a pain for something like
HMMs because the unobserved discrete states need to be marginalized
out.

- Bob

> On Jul 17, 2015, at 1:42 AM, Wad Wan <cphni...@gmail.com> wrote:
>
>
> Yo-yo, check my rhyme!..

>
> So, guys to be more specific, I mean the following. Suppose that we have a financial or economic time series, where we have a prior subjective belief that
> there are two different latent states (i.e., bull and bear markets). Now, we want to use a hidden Markov models (i.e., Markov Regime State Switching) that is
> somewhat similar to SWARCH model by Hamilton and Susmel. But to keep it simple and stupid, we begin to look at ARCH(1) model and Markov Regime State Switching Model. Now, what I have a great difficulty to figure out how to use Stan to combine these two models together. When I read Stan's user manual, they
> write very well the explanation for how to solve ARCH(1) model, and the explanation for how to solve HMM, but there is that missing link that how do you combine these two models in Stan, such that, we run HMM recursively through the time series, and use ARCH(1) model to optimize parameters in the two latent states?...
>
> So, to keep it simple and stupid we begin with the first order ARCH(1) stochastic process estimation of paramters in the model. The definition of ARCH(1) looks something like that:
>
> \begin{eqnarray}
> a_t = \sigma_t \epsilon_t, \text{ } \sigma_t = \sqrt[]{\alpha_0 + \alpha_1 a^2_{t-1}}, \text{ } \sigma_t^2 = \alpha_0 + \alpha_1 a^2_{t-1},
> \end{eqnarray}
>
> where $\alpha_0 > 0$, $\alpha_1 \geq 0$, $\epsilon_t$ is independent of the past, and $\sigma_t$ is known function value of the prior information
> at time $t-1$.
>
> We begin our implementation in R and Stan, which looks like that:
>

> # 1st we simulate ARCH(1) stochastic process, with known user-defined-true-parameters
>

> library(rstan) # we load in the Stan package
> library(tseries)
>

> # In first step, we generate ARCH(1) process
> #

> alpha0 <- 0.06321 # user defined subjective constant
> alpha1 <- 0.03548 # user defined subjective constant
> T <- 10000 # number of replications
>
> r_t <- rep(0,T)
>
> sigma2_t <- 0.1 # initial guess; user defined parameter
> r_t[1] <- sqrt(sigma2_t)*rnorm(1,0,1)
>
> for (i in 2:T) {
> eps <- rnorm(1,0,1)
> sigma2_t <- alpha0 + alpha1*r_t[i-1]^2
> r_t[i] <- sqrt(sigma2_t)*eps
> }
>
> plot.ts(r_t)

> hist(r_t)

>
> ## Stan modelling in R
> #

> modelString = "
>
> data {
> int<lower=0> T; // number of time points
> real r_t[T]; // return at time t
> }
>
> parameters {
> real mu;
> real<lower=0> alpha0; // noise intercept
> real<lower=0,upper=1> alpha1; // noise slope
> }
>
> model {
>
> for (t in 2:T)
> r_t[t] ~ normal(mu, sqrt(alpha0 + alpha1 * pow(r_t[t-1] - mu,2)));
> }
> "
>
> sm <- stan_model(model_code = modelString) # Stan model
>
> fit <- sampling(sm, data=list(r_t=r_t, T=T), iter=1000, chains=4, warmup=250, thin=5); fit
>
> print(fit, digits_summary = 5)
> plot(fit)
>
> fit_sim <- extract(fit)
> hist(fit_sim$alpha0)
> hist(fit_sim$alpha1)

> # Done!..

>
> Den onsdag den 15. juli 2015 kl. 01.00.59 UTC+2 skrev Wad Wan:
> Hello guys,
>
> I want to do ARCH(1) model for 2-state Markov Regime Switching Model for some financial time series.
>
> In Stan Modeling Language User's Guide and Reference Manual (v.2.6.2) on page 69
> I have found how to implement ARCH(1) model. This code works well, and I get correct approximations.
>
> On page 81 I find the HMM "Semisupervised Estimation", but I can't figure out how to combine the ARCH(1)
> on page 61 and Semisupervised Estimation on p. 81.
>
> Does anyone have an idea how to do it?..
>
> All the best, Serge
>

> --
> 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/d/optout.

[andre.p...@googlemail.com]

Serge,

suppose we at the time t, then we would need to update the calculation of the ARCH parameters for all states. 

With updating the sequence of states may chance. Thus in every iteration we need to calculate the path of

sequences of the states. However this would result in discontinuous objective function.

Just my 2 cents,

Andre

[Wad Wan]

Hey Andre,

Thank a lot for your 2 cents, it is so cool that want to give some feedback about it. 

So, check it out, it goes something like this. I have found the solution for ARCH(1) with jumps in regime, that is already solved by Hamilton and Susmel (1994). But the thing, is it is implemented in GAUSS programming language. I am trying to translate the program from GAUSS to R, but it is not that easy as I thought it was. Now, soon I am going to do debugging in R, but I don't really know how to. Because, I have two functions, and I need to have one main function that calls these two functions. Do you, by any chance, know how to do debugging of the main function in R (if there are any)?

[Wad Wan]

Den onsdag den 15. juli 2015 kl. 01.00.59 UTC+2 skrev Wad Wan:

By the way, the guys from Stan-team, haven't implemented the autoregressive processes with discrete jumps, that is why it is impossible to make it work in Stan. I am working currently at the translation of Hamilton and Susmel (1994) model, but the main challenge is the debugging of it. Does anybody can help me, cause I'm kinnda stop in the middle. The thing is that R's vector definition and operations is a bit different from GAUSS programming language.  

[Bob Carpenter]

I think such a model could be implemented in Stan if I understand
you correctly --- there are multiple states, and the observation
is a function of both the current state (like any state-space model)
and in this case, some history. Not easy to code, though.

- Bob

> On Oct 18, 2015, at 2:02 AM, Wad Wan <cphni...@gmail.com> wrote:
>
>
>
> Den onsdag den 15. juli 2015 kl. 01.00.59 UTC+2 skrev Wad Wan:
> Hello guys,
>
> I want to do ARCH(1) model for 2-state Markov Regime Switching Model for some financial time series.
>
> In Stan Modeling Language User's Guide and Reference Manual (v.2.6.2) on page 69
> I have found how to implement ARCH(1) model. This code works well, and I get correct approximations.
>
> On page 81 I find the HMM "Semisupervised Estimation", but I can't figure out how to combine the ARCH(1)
> on page 61 and Semisupervised Estimation on p. 81.
>
> Does anyone have an idea how to do it?..
>
> All the best, Serge
>
>

> By the way, the guys from Stan-team, haven't implemented the autoregressive processes with discrete jumps, that is why it is impossible to make it work. I am working currently at the translation of Hamilton and Susmel (1994) model, but the main challenge is the debugging of it. Does anybody can help me, cause I'm kinnda stop in the middle. The thing is that R's vector definition and operations is a bit different from GAUSS programming language.

>
> --
> 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.

> To post to this group, send email to stan-...@googlegroups.com.

[ibby]

Hello Bob,

The idea is presented in: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.453.8556&rep=rep1&type=pdf

 

If there is code for HMM, it should be easy to implement this system.

[Bob Carpenter]

Yup --- that's about what I would have expected. And yes, you
should be able to implement it in Stan.

- Bob