I want to draw Bayesian inference from a time-series of univariate daily logreturns, denoted as $y\in\mathbb{R}^T$. My likelihood is chosen to be t-distributed with unknown parameters $\mu,\sigma^2,\nu$: $$\log\mathcal{L}(y|\mu,\sigma^2,\nu)=T[\log\left(\frac{\Gamma(\nu+1)/2}{\Gamma(\nu/2)}\right)-0.5\log(\nu\pi)-0.5log(\sigma^2)]\\-(\nu+1)/2\sum\limits_{i=1}^{T}log\left(1+\frac{(y-\mu)}{\nu\sigma^2}^2\right).$$
I do not have subjective beliefs regarding any of the parameters. Therefore I am thinking of an flat prior $p(\mu)\propto c$, Jeffreys prior $p(\sigma^2)\propto \frac{1}{\sigma^2}$ and $p(\nu)\propto \frac{1}{\nu}$. I want to compare the effect of choosing importance densities for $\nu$, namely: $$p_1(\nu)\sim\mathcal{U}[0,\nu_\max]\\p_2(\nu)\sim\mathcal{E}(\lambda).$$ where $\mathcal{E}$ denotes the exponential prior. The choice of $p_1(\nu)$ comes from this post, whereas $p_2(\nu)$ and $p_3(\nu)$ are standard ideas. $\nu_\max$ is chosen to be 50, as above a certain threshold I'd simply conclude that the data follows a normal distribution. I set $a$ to be the MLE estimator and $\lambda$ is $4$, a choice that also reflects some assumptions regarding the true parameter.
I try to implement these different models via an importance sampling algorithm. However, I think that my implementation contains either an error or my choice of the importance density is in some way not sensitive enough. In the following code I generate artificial data from a $t$ distribution with known parameters and compute bayesian estimated for each of the mentioned importance densities. Based on the created sample I repeated the same analysis $100$ times to obtain some rough idea of the variance of my estimates. The results for the degree of freedom are illustrated in this boxplot (the horizontal line shows the MLE estimator for $\nu$)
I am surprised by the huge variance implied by the second importance density compared to the first one. However, as there are quite a few parameters I had to set arbitrary (number of importance density draws, $a$, $\lambda$, etc) I would like to know how to 'robustify' my results: Are there any measures to check if my posterior mean estimates of $\nu$ are appropriate? I appreciate each and every comment on how to use the available data to argue in favor of one or another model.
Update: I just figured out that the resulting weights in my application are in roughly $99$ percent of the cases equal to 0. This certainly comes from a poor choice of my importance densities. Based on my setting, is this obvious a priori?
close all;clc;
s = rng; %Set seed
T=500;
data=mvtrnd(1,4,T);
logL_t=@(y,para) length(y)*(log(gamma((para(3)+1)./2))-log(gamma(para(3)...
/2))-0.5*log(para(3)*pi)-0.5*log(para(2)))-(para(3)+1)/2....
*sum(log(1+1./(para(3)*para(2)).*(y-para(1)).^2));
logL_n=@(y,para) -length(y)/2*log(2*pi*para(2))-1/(2*para(2))*...
sum((y- ara(1)).^2);
logL_e=@(y,para) log(para)-para*y;
logL_sigma=@(y,para) log(1/(2*para)*(y<para)+1/(2*para)*(y>para)*y^(-2));
[T,~]=size(data);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% ML Properties %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t_mu_0 = 8;t_var_0 = 3;df_0 = 20;
Init_point = [t_mu_0 t_var_0 df_0]';
[MLE,~,~,~,~,hessian]=fminunc(@(val_param) -1*logL_t(data,val_param),Init_point);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Importance Sampling %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j=1:100
display(j)
N_S = 10000; %Number of proposals
prop_var = 5; %Variance of the parameter mu
prop_mean=0; %Mean of the parameter mu
prop_lambda=4; %Exponential distribution parameter for nu_3
U_max=30; %Parameter for proposal of nu_2
posterior1= zeros(N_S,1);posterior2= zeros(N_S,1);posterior3= zeros(N_S,1); %Define empy matrices
prop_dens1= zeros(N_S,1);prop_dens2= zeros(N_S,1);prop_dens3= zeros(N_S,1);
log_likeli1 = zeros(N_S,1);log_likeli2 = zeros(N_S,1);log_likeli3 = zeros(N_S,1);
prop_d_mu=normrnd(prop_mean,sqrt(prop_var),N_S,1); %Proposition for mu
prop_d_sigma=0.5*rand(N_S,1)*sqrt(MLE(2))+0.5* (rand(N_S,1)*1/sqrt(MLE(2))).^(-1); %Proposition for sigma
prop_d_v2=rand(N_S,1)*U_max; %Proposition for nu_2
prop_d_v3=exprnd(prop_lambda,N_S,1); %Proposition for nu_3
prop_draw2=[prop_d_mu,prop_d_sigma,prop_d_v2];
prop_draw3=[prop_d_mu,prop_d_sigma,prop_d_v3];
for i=1:N_S
param2 = prop_draw2(i,:);
param3 = prop_draw3(i,:);
log_likeli2(i) = logL_t(data,param2);
log_likeli3(i) = logL_t(data,param3);
prop_dens2(i)= logL_n(param2(1),[prop_mean,prop_var])- logL_sigma(param2(2),MLE(2))-log(U_max);
prop_dens3(i)= logL_n(param3(1),[prop_mean,prop_var])- logL_sigma(param3(2),MLE(2))+logL_e(param3(3),prop_lambda);
posterior2(i) = log_likeli2(i)-log(param2(2))-log(param2(3));
posterior3(i) = log_likeli3(i)-log(param3(2))-log(param2(3));
end
diff_max2 = max(posterior2-prop_dens2);
diff_weight2 = exp(posterior2-diff_max2);
log_diff_cst2 = log(mean(diff_weight2)) + diff_max2; % Normalizing constant
diff_weight2 = diff_weight2/sum(diff_weight2); % Weight for each particle
diff_max3 = max(posterior3-prop_dens3);
diff_weight3 = exp(posterior3-diff_max3);
log_diff_cst3 = log(mean(diff_weight3)) + diff_max3; % Normalizing constant
diff_weight3 = diff_weight3/sum(diff_weight3); % Weight for each particle
r(j,:,1)=sum(prop_draw2.*repmat(diff_weight2,1,3));
r(j,:,2)=sum(prop_draw3.*repmat(diff_weight3,1,3));
end
boxplot(reshape(r(:,3,:),j,2))
hold on
plot(xlim,[MLE(3),MLE(3)])
title('Bayesian estimates after 100 repetitions')
