Hi - I’m fitting a model for an error structure for a structural vector autoregression, but I the issue I’m encountering is much more basic.
Essentially, for data y, my model has
y[,t]=H*diag_matrix(vol[,t])*z[,t],
that is, y is a mixture of white noise series z with stochastic volatilities. I can implement this easily enough by using an ugly loop with many inversions in my model block:
for (t in 1:T) y[,t]~multi_normal(zeros,H*diag_matrix(vol[,t].*vol[,t])*H')
This runs slowly. It seems to me that a more desirable approach would be to express (presumably in transformed parameters block?)
eps=inv(H)*y zhat=eps./vol
and then, in my model block, I have
zhat[i,]~normal(0,1) for each dimension i. Unsurprisingly, this runs much faster (like two orders of magnitude faster), as I’m avoiding both a loop and the multi_normal. However, the results aren’t “right”. By that I mean, for simulated data the first approach estimates the true parameters well, the latter does terribly. Mathematically, I’m doing the same thing - just standardizing a random normal before computing the likelihood. The rest of the code is identical, I just change the element whose value enters the likelihood in the model block. Any ideas what’s causing the dramatic difference in accuracy? Comparing to other things I’ve run, it looks to me that Stan is possibly ignoring the distribution of zhat specified in the model completely. I’m attaching sample code, in case that makes things clearer (or you catch typos!). Thanks so much!
Daniel
SV_code= {
'data {'
' int<lower=0> T; // # time points (equally spaced)'
'matrix[2,T] y;
'}'
'transformed data {'
'vector[2] zeros; '
'zeros = rep_vector(0, 2);'
'}'
'parameters {'
'vector<lower=-1,upper=1>[2] phi;'
' cov_matrix[2] Sigma; // innovation variance'
' matrix[2,T] h_innov;// std log variance time t'
'real H21;'
' real H12;'
'}'
'transformed parameters {'
' matrix[2,T] h; // log variance at time t'
'matrix[2,2] H;'
'matrix[2,2] L;'
'matrix[2,T] zhat;//comment out for other method'
'L=cholesky_decompose(Sigma);'
'H[1, 1] = 1; H[1, 2] = H12; // Construct normalized H matrix'
'H[2, 1] = H21; H[2, 2] = 1;'
' h = L*h_innov; // Assemble the path of volatility'
' h[1,1] = h[1,1] / sqrt(1 - phi[1] * phi[1]);'
' h[1,2] = h[1,2] / sqrt(1 - phi[2] * phi[2]);'
' for (t in 2:T){'
' h[,t] = h[,t] +diag_matrix(phi) * h[,t-1];'
'}'
'zhat=inv(H)*y; // standardize the data,y to something that is iid N(0,1) //comment out for other method'
'zhat=zhat ./ exp(h/2);//comment out for other method'
'}'
'model {'
'zhat[1,]~normal(0,1);//comment out for other method'
'zhat[2,]~normal(0,1);//comment out for other method'
'h_innov[1,]~normal(0,1);'
'h_innov[2,]~normal(0,1);'
'for (t in 1:T){'
'//y[,t]~multi_normal(zeros,H*diag_matrix(exp(h[,t]))*H'');// uncomment this line to run without standardization'
'}'
'}'
};
rng(1)
T=500;
H=[1,-.5;.5,1];
innov=[0.2,-0.025;-0.025,0.2];
innovchol=chol(innov);
eta=NaN(2,T);
X=sqrt(diag(innov))./sqrt(1-phi.^2).*randn(2,1);
V=exp(X);
z=sqrt(V).*randn([2 1]);
eta(:,1)=H*z;
for t=2:T
X=phi.*X+innovchol'*randn([2,1]);
V=exp(X);
z=sqrt(V).*randn([2 1]);
eta(:,t)=H*z;
end
SV_data = struct('T',length(eta),'y',eta);
fitcomplile=stan('model_code',SV_code,'data',SV_data,'warmup',10,'iter',10,'chains',1);
fit=stan('fit',fitcompile,'data',SV_data,'warmup',100,'iter',1000,'chains',1,'seed',1);
%%% With standardization, H estimated to be [1.0000 12.2478 ;-17.4181
%%% 1.0000]. With the multi_normal, [1.0000 1.3922; -1.9250
%%% 1.0000] which is close to a column-reordering of H given the small
%%% run-length