Hi all,
I am new to Stan and have trouble estimating process and observation error standard deviations using an extended Kalman filter (EKF). Following Durbin & Koopman (2012), one can use Maximum likelihood estimation to estimate parameters in the process and observation equations. In a frequentist version, this works reasonably well. Now, I wanted to implement a Bayesian version. In this simplest example possible, I only try to estimate the two standard deviations, but the model notoriously overestimates both substantially (should be 0.05 both, but estimates are far above). Here, I use a normal prior, but I also tried a lognormal prior with similar issues. The current priors are pretty uninformative, but it also happens with really informative priors.
I would like some advice as to why this problem persists, or whether I made a coding error that causes this?
The (simulated) data stems from the R code here below, I use rstan to call the Stan code named “BESAstan_onlysig.stan”, which is attached here below.
All the best and thank you very much in advance,
Kira
library(rstan)
# To run in parallel on multiple cores
options(mc.cores=parallel::detectCores())
# To avoid recompiling unchanged Stan program
rstan_options(auto_write=TRUE)
## parameters, then simulate data (pt and Ht time series, time series of length N)
a=0.5737
b=0.87
c=2
d=0.7
r=0.3809
K=500
kappa=r/K
sigmaeps=0.05;
sigmaeta=0.05
sigmaphi=0.1
N=50;
psi=log(250)
xt<- vector(length=N)
xt[1]=psi
epsilont=0
etat=0
phit=0
pt=0
ht=0
for (j in 1:N) {
epsilont[j]=rnorm(1,0,sigmaeps)
etat[j]=rnorm(1,0,sigmaeta)
phit[j]=rnorm(1,0,sigmaphi)
pt[j]=log(c)-d*xt[j]+epsilont[j];
ht[j]=1/b*(log(a)-pt[j]+phit[j]);
xt[j+1]=log((1+r)*exp(xt[j])-kappa*exp(xt[j])^2-exp(ht[j]))+etat[j];
}
Ht=exp(ht) # needed in levels in the filter
# define prior distribution parameters:
musigeps=1
sigsigeps=1
musigeta=1
sigsigeta=1
nsamples=6000
nchains=3
# feed in data and run model
dataStan <- list(pt=pt,Ht=Ht, TT=N,mur=r,musigeps=musigeps,sigsigeps=sigsigeps,musigeta=musigeta,sigsigeta=sigsigeta)
f2pStan <- stan(file = "BESAstan_onlysig.stan", data = dataStan, chains = nchains, iter =nsamples,warmup=1000,thin=2)
# Traceplot (includes warm up )
traceplot(f2pStan, pars="sigeps", inc_warmup = TRUE)
traceplot(f2pStan, pars="sigeta", inc_warmup = TRUE)
# print some summary results: true means of sigmaeps and sigmaeta should be around 0.05!
print(f2pStan, max = 150)
/*----------------------- Data --------------------------*/
/* Data block: unpack the objects that will be inputted as data */
data {
int TT; // Length of state and observation time series
vector[TT] pt; // observations (left side of observation equation)
vector[TT] Ht; // observations (used in process equation)
real mur; // parameter in process equation
real musigeps;
real sigsigeps;
real musigeta;
real sigsigeta;
}
/*----------------------- Parameters --------------------------*/
parameters {
real<lower=0> sigeps; // Standard deviation of the observation equation
real<lower=0> sigeta; // Standard deviation of the process equation
}
/*----------------------- Model --------------------------*/
/* Model block: defines the model */
model {
// Initialization of variables
vector[TT+1] st; // state variable
vector[TT] vt; // residual time series
vector[TT] Ft; // variance of vt
vector[TT+1] Zt; // variance of st
vector[TT] stt; // intermed state var time series
vector[TT] Ztt; // intermed state var variance
vector[TT] Jtt; // Jacobian
real llik_obs[TT];
real llik;
// Prior distributions:
target+=normal_lpdf(sigeps|(musigeps),sigsigeps);
target+=normal_lpdf(sigeta|(musigeta),sigsigeta);
// initial:
Zt[1]=0; // initial variance
st[1]=log(0.5); // initial state
target+= -0.5*TT*log(2*3.14159); // log likelihood constant part
// Extended Kalman filter part loop
for (t in 1:TT){
vt[t]=pt[t]-log(2)+0.7*log(500)+0.7*st[t]; // residual
Ft[t]= pow(0.7,2)*Zt[t]+pow(sigeps,2); // variance of residual
// Each step, add this part to loglikelihood (Durbin & Koopmann 2012, p. 171):
llik_obs[t] = -0.5*(log(Ft[t]+pow(vt[t],2)/Ft[t]));
stt[t]=st[t]-(0.7*Zt[t]*vt[t])/Ft[t]; // corrected state estimate
Ztt[t]=pow(sigeps,2)*Zt[t]/Ft[t]; // variance of corrected state estimate
st[t+1]=log((1+mur)*exp(stt[t])-mur*exp(2*stt[t])-Ht[t]/(500)); // a priori state estimate for next period
Jtt[t]=((1+mur)*exp(stt[t])-2*mur*exp(2*stt[t]))/((1+mur)*exp(stt[t])-mur*exp(2*stt[t])-Ht[t]/(500)); // Jacobian
Zt[t+1]=pow(Jtt[t],2)*Ztt[t]+pow(sigeta,2); // a priori estimate for state variance for next period
}
llik = sum(llik_obs); // log likelihood is the sum over all the parts
target+=llik; // add the log likelihood to the target
}