function [S,M,V]=log_chisquared1_approx_states(Z)

%  Compute the matrix containing the indexes of components of the mixture
%  of normal distributions approximating the log_chi_squared_1 
%  as in Kim, Shephard, Chib (1998)
%
%  If   z_t is a mx1 vector of independent log_chi_squared_1 distributions,
%  and z_t is independent of z_(t+j) for all j, then given a the
%  realizations of z_t we can draw the specific components from the
%  following discrete distribution
%  P(s_t(h)=j)= m_p(j) * normpdf( z_t(h),m_m(j), m_v(j) )  
%  for all j=1:7, for all h=1:m and all t=1:T
%  where P has to be normalized (normalization will be done by a subroutine)
%
%  Update (October 2014): the mixture is updated to 10 components, 
%  following Omori, Chib, Shephard and Nakajima (2007). 
%
%  INPUTS: Z is mxT
%  OUTPUTS:  S is mxT, contains the drawn indexes of components
%            M is mxT, contains the associated means of components
%            V is mxT, contains the associated variances of components
%           

% Copyright (c) 2013-2018, Francesco Corsello, All Rights Reserved.

[m,T]=size(Z);


% OLD MIXTURE 
% weights of the mixture
%m_p=[7.3e-3,.10556,2e-5,.04395,.34001,.24566,.2575];
%moments of the normal components
% -1.2704  is added to the mean vector to match the first moment  
%m_m=[-10.12999,-3.97281,-8.56686,2.77786,0.61942,1.79518,-1.08819]-1.2704;
%m_v=[5.79596,2.61369,5.17950,0.16735,0.64009,0.34023,1.26261];


%NEW MIXTURE
m_p=[ 0.00609, 0.04775, 0.13057, 0.20674, 0.22715, 0.18842, 0.12047,...
    0.05591, 0.01575, 0.00115];
m_m=[ 1.92677, 1.34744, 0.73504, 0.02266, -0.85173, -1.97278, -3.46788,...
    -5.55246, -8.68384, -14.65000];
m_v=[ 0.11265, 0.17788, 0.26768, 0.40611, 0.62699, 0.98583, 1.57469,...
     2.54498, 4.16591, 7.33342];
 

%stack in matrices
vZ=repmat(Z(:),1,length(m_p));
v_weights=repmat(m_p,m*T,1);
v_mean=repmat(m_m,m*T,1);
v_sigma=sqrt(repmat(m_v,m*T,1));

% compute the probabilities
P=v_weights .* normpdf( vZ ,v_mean, v_sigma );
Svec = sample_discr_same_length(P);
S=reshape(Svec,m,T);

m_vec=m_m(Svec);
M=reshape(m_vec,m,T);
var_vec=m_v(Svec);
VV=reshape(var_vec,m,T);
V=cell2mat(reshape(arrayfun(@(t) diag(VV(:,t)),1:T,'UniformOutput',false),1,1,T));
  



