% MCMC experiment to simulate the VIX return distribution

load('VIX_return_data');  % load the data set for returns
nstate = 10;              % number of states

% remove some outliers:
vixret = vixret( vixret < 0.4 );
vixret = vixret( vixret > -0.2);

% target distribution to be simulated:
pai = hist(vixret,nstate)/length(vixret);

plot(vixret(length(vixret):-1:1),'*');
xlabel('time'); ylabel('daily return');

Q=zeros(nstate,nstate);
bin_size = 0.75;          % "optimal" bin size based on experiments

% compute the transition probability matrix, notice that it's not
% symmetric.
for i=1:nstate
    for j=1:nstate
        if j==1
            yr = normcdf((j-i+0.5)*bin_size);
            yl = 0;
        elseif j==nstate
            yr = 1;
            yl = normcdf((j-i-0.5)*bin_size);
        else
            yr = normcdf((j-i+0.5)*bin_size);
            yl = normcdf((j-i-0.5)*bin_size);
        end
        Q(i,j) = yr - yl;
    end
end

N_sim = 5000;          % number of realizations
N_step = 100;          % number of steps in the Markov chain
end_state = zeros(N_sim,1);

for k=1:N_sim
    state = 8;         % start from state "4"
    for j=1:N_step
    state_next = update_state(state,bin_size);
    alpha = min(pai(state_next)*Q(state,state_next)/pai(state)/Q(state_next,state),1);
    if rand() < alpha
        state = state_next;
    end
    end
    end_state(k) = state;
end
hist(end_state,10);
title('simulated end states histogram');
    
    

