Hello,
I am trying to fit a model known as Autoregressive Conditional Duration (ACD) in Stan. ACD model is very similar to GARCH.
x_i=\psi_i\epsilon_i
\psi_i=\omega+\alpha x_{i-1}+\beta \psi_{i-1}
\epsilon_i \sim exponential(1)
To verify the quality of fit I simulate from a known model first and then at a later will use real data.
Attached is the code to simulated data and stan model code
#ACD- 1 regime
x<-c()
psi<-c()
epsilon<-c()
N<-2000
epsilon<-rexp(N,1)
omega<-0.54
alpha<-0.45
beta<-0.25
psi[1]<-omega/(1-alpha-beta)
x[1]<-psi[1]*epsilon[1]
for(i in 2:N)
{
psi[i]<-omega+alpha*x[i-1]+beta*psi[i-1]
x[i]<-psi[i]*epsilon[i]
}
data<-list(N=length(x),x=x)
modelcode_acd<-"
data{
int<lower=0> N;
vector[N] x;
}
parameters{
real<lower=0> omega;
real<lower=0, upper=1> alpha;
real<lower=0, upper=(1-alpha)> beta;
}
transformed parameters{
real<lower=0> psi[N];
psi[1]=omega/(1-alpha-beta);
for(i in 2:N)
psi[i]=omega+alpha*x[i-1]+beta*psi[i-1];
}
model{
real one_over_psi;
for(i in 1:N){
one_over_psi=1/psi[i];
target+=exponential_lpdf(x[i] | one_over_psi);
}
}"
acd_fit <- stan(model_code = modelcode_acd, data = data)
I am getting results which are fair enough. True values of (\omega,\alpha,\beta) is (0.54,0.45,0.25). Here is the summary of result that I am getting:
Inference for Stan model: f94ed5db5c51752f9e252f9e086fc4f1.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean | se_mean | sd | 2.50% | 25% | 50% | 75% | 97.50% | n_eff | Rhat | |
---|---|---|---|---|---|---|---|---|---|---|
alpha | 0.49 | 0 | 0.03 | 0.42 | 0.47 | 0.49 | 0.51 | 0.55 | 2263 | 1 |
beta | 0.17 | 0 | 0.04 | 0.09 | 0.15 | 0.17 | 0.2 | 0.26 | 1470 | 1.01 |
omega | 0.62 | 0 | 0.06 | 0.51 | 0.58 | 0.62 | 0.66 | 0.74 | 1430 | 1.01 |
My questions for the above is: Is there any way I can further improve the results?
The above model was a one regime model. I want to further try a two regime model in which the change is as follows:
\psi_i=\omega_i+\alpha_i x_{i-1}+\beta_i \psi_{i-1}
(\omega_i,\alpha_i,\beta_i)= (\omega_k,\alpha_k,\beta_k)~ with ~probability~ p_k~\forall i
k=1~ and~ 2~ with ~p_1+p_2=1
Attached is the simulation code and model code.
#ACD-2 regime
N<-2000
x<-c()
psi<-c()
epsilon<-c()
epsilon<-rexp(N,1)
omega<-c(0.54,0.3)
alpha<-c(0.35,0.55)
beta<-c(0.15,0.25)
pi<-c(0.3,0.7)
z <- vector("numeric", N)
z<- sample(1:2, size = N, replace=TRUE, prob = pi)
psi[1]<-omega[z[1]]/(1-alpha[z[1]]-beta[z[1]])
x[1]<-psi[1]*epsilon[1]
for(i in 2:N)
{
psi[i]<-omega[z[i]]+alpha[z[i]]*x[i-1]+beta[z[i]]*psi[i-1]
x[i]<-psi[i]*epsilon[i]
}
data<-list(N=length(x),x=x)
modelcode_acd<-
"data{
int<lower=0> N;
vector[N] x;
}
parameters{
simplex[2] pi;
positive_ordered[2] omega;//for identifiability
real<lower=0, upper=1> alpha[2];
real<lower=0, upper=1-alpha[1]> beta1;
real<lower=0, upper=1-alpha[2]> beta2;
}
transformed parameters{
real<lower=0> beta[2];
vector<lower=0>[2] psi[N];
beta[1]=beta1;
beta[2]=beta2;
psi[1,1]=omega[1]/(1-alpha[1]-beta[1]);
psi[1,2]=omega[2]/(1-alpha[2]-beta[2]);
for(i in 2:N){
for(k in 1:2){
psi[i,k]=omega[k]+alpha[k]*x[i-1]+beta[k]*psi[i-1,k];
}
}
}
model{
vector[2] contributions;
real one_over_psi;
pi ~ dirichlet(rep_vector(1.0,2));
// likelihood
for(i in 1:N) {
for(k in 1:2) {
one_over_psi=1/psi[i,k];
contributions[k] = log(pi[k]) + exponential_lpdf(x[i] | one_over_psi);
}
target += log_sum_exp(contributions);
}
}"
acd_fit <- stan(model_code = modelcode_acd, data = data, chains=4, iter=10000,control=list(adapt_delta=0.9))
For now, I am not interested in estimating regime for each observation and would like to restrict to estimating predictive density only. But, if at all, that is possible it will be definitely an advantage.
Running the above code gives a lot of divergent transition. There is no point in accessing the quality of the fit. Is there any way to decrease the number of divergent transitions? Or shall I ignore this warning and continue my further modeling?
Once I am able to get over this problem, I want to generalize this from 2 to K regimes where K is any finite number.
I have read relevant blogs(or articles) on mixture modeling, Stan and mixture modeling in Stan.
Kindly suggest.