 function [logL,likelihoods ,rho, smoothed]=P_RS_dcc_likelihood(param, data,isnorm)
%说明
%输入变量：
%para 模型参数。顺序[a1 b1 g1 a2 b2 g2 p11 p22 nu1 nu2],符号含义与模型对应；
%data 数据。结构[stdresid，ht，D]，stdresid 为garch 拟合后的标准化残差，ht
%是garch拟合后的标准差，D是非对称性冲击变量；
%isnorm 取 1或0，1表示为正态分布，0表示二元t分布；
%输出变量：
%logL 对数似然值
%likelihood 似然值
%rho 相关系数；
%smoothed 平滑概率


stdresid = data(:,1:2);%标准化后的残差
ht = data(:,3:4);%标准差
%D = data(:,4:5);%非对称性变量


[T,k] = size(stdresid);
Qbar = cov(stdresid); %%初始值设定，此处为状态1，表示经济形式好的时候
Qbar2 = cov(stdresid)+0.05*(ones(k,k)-[1 0;0 1]);%状态2，表示经济形势差，相关性增大
%Nbar = cov(D);%状态1
%Nbar2 = cov(D)-0.03*(ones(k,k)-[1 0;0 1]);%状态2

Q00 = zeros(k,k,T);%00表示状态1
R00 = zeros(k,k,T);%R=sqrt(diag(Q))^(-1)*Q*sqrt(diag(Q))^(-1)
H00 = zeros(k,k,T);%H=h_t*R_t*h_t
d00 = zeros(T,1);%密度函数值

Q11 = zeros(k,k,T);%11表示状态2
R11 = zeros(k,k,T);
d11 = zeros(T,1);
H11 = zeros(k,k,T);

Q00(:,:,1) = Qbar;%令Q的初始值为Qbar
Q11(:,:,1) = Qbar2;

p00 = zeros(T,1);%对应hamilton 书中的 xi_t|t_1
p11 = zeros(T,1);
p0 = zeros(T,1);%xi_t|t
p1 = zeros(T,1);

a1 = param(1);
b1 = param(2);
%g1 = para(3);
a2 = param(3);
b2 = param(4);
%g2 = para(6);
transP11 = exp(param(5))/(1+exp(param(5)));%由状态1转到状态1的概率
transP22 = exp(param(6))/(1+exp(param(6)));%由状态2转到状态2的概率

%判断是否为t分布，取t分布的参数值
if isnorm == 0
    nu1 = param(7);
    nu2 = param(8);
end

P = [transP11 1-transP11;1-transP22 transP22];%形成一个转移概率矩阵

p0(1,1) = (1-transP22)/(2-transP11-transP22);%无条件的概率初始值设定，参考hamilton设定方式
p1(1,1) = (1-transP11)/(2-transP11-transP22); 

p00(1:2,1) = [p0(1,1);p0(1,1)];
p11(1:2,1) = [p1(1,1);p1(1,1)];

rho = zeros(T,1);%相关系数


logL=0;%对数似然值
likelihoods=zeros(T,1);

stdresid=[zeros(1,k);stdresid];
for j=2:T
    %% 状态1 的密度密度函数值计算
   Q00(:,:,j) = Qbar*(1-a1-b1);%-Nbar*g1; %常数部分  
   Q00(:,:,j) = Q00(:,:,j)+a1*(stdresid(j-1,:)'*stdresid(j-1,:));%加上残差部分
   Q00(:,:,j) = Q00(:,:,j)+b1*Q00(:,:,j-1);%加上滞后一项
%   Q00(:,:,j) = Q00(:,:,j)+g1*(D(j-1,:)'*D(j-1,:));%加上非对称性冲击
   R00(:,:,j) = Q00(:,:,j)./(sqrt(diag(Q00(:,:,j)))*sqrt(diag(Q00(:,:,j)))');
   H00(:,:,j) = diag(ht(j,:))*R00(:,:,j)*diag(ht(j,:)); %Ht=ht*Rt*ht
   if isnorm == 0 %计算二元t分布密度函数值
      d00(j) = gamma((nu1+k)/2)./(gamma(nu1/2)*(pi*(nu1))^(k/2)*sqrt(det(H00(:,:,j))))*(1+stdresid(j,:)*inv(R00(:,:,j))*stdresid(j,:)'/(nu1))^(-nu1/2-k/2);
   else %计算二元正态分布密度函数值
      d00(j) = 1./sqrt(det(H00(:,:,j))*(2*pi)^(k))*exp(-1/2*stdresid(j,:)*inv(R00(:,:,j))*stdresid(j,:)');
   end
   %%
   %状态2的密度函数值计算，过程和状态1一致
   Q11(:,:,j) = Qbar2*(1-a2-b2);%-Nbar2*g2;   
   Q11(:,:,j) = Q11(:,:,j)+a2*(stdresid(j-1,:)'*stdresid(j-1,:));
   Q11(:,:,j) = Q11(:,:,j)+b2*Q11(:,:,j-1);
%   Q11(:,:,j) = Q11(:,:,j)+g2*(D(j-1,:)'*D(j-1,:));
   R11(:,:,j) = Q11(:,:,j)./(sqrt(diag(Q11(:,:,j)))*sqrt(diag(Q11(:,:,j)))');
   H11(:,:,j) = diag(ht(j,:))*R11(:,:,j)*diag(ht(j,:)); %Ht=ht*Rt*ht
   if isnorm == 0
      d11(j) = gamma((nu2+k)/2)./(gamma(nu2/2)*(pi*(nu2))^(k/2)*sqrt(det(H11(:,:,j))))*(1+stdresid(j,:)*inv(R11(:,:,j))*stdresid(j,:)'/(nu2))^(-nu2/2-k/2);
   else 
      d11(j) = 1./sqrt(det(H11(:,:,j))*(2*pi)^(k))*exp(-1/2*stdresid(j,:)*inv(R11(:,:,j))*stdresid(j,:)');
   end
   %%
   likelihoods(j) = d00(j,1)*p00(j,1)+d11(j,1)*p11(j,1); %% 总的似然值=状态1似然值乘以对应的概率+状态2似然值乘以概率 
   logL=logL+log(likelihoods(j));%对数似然值，并加总
   
   rho(j,1) = R00(1,2,j)*p00(j,1)+R11(1,2,j)*p11(j,1); %%总的相关系数，计算方法如总的似然值

   p0(j+1,1) = d00(j,1)*p00(j,1)/likelihoods(j); %更新概率
   p1(j+1,1) = d11(j,1)*p11(j,1)/likelihoods(j); %

   p00(j+1,1) = P(1,1)*p0(j+1,1)+P(1,2)*p1(j+1,1);% 用转移概率更新滤波概率，即xi_tt_1=P*xi_tt
   p11(j+1,1) = P(2,1)*p0(j+1,1)+P(2,2)*p1(j+1,1);
end;

rho = rho(2:end);

logL = -logL;
likelihoods = likelihoods(2:T);

%%
%以下计算平滑概率
s00 = zeros(T,1);
s11 = zeros(T,1);
s00(T,1) = p0(T,1);%赋于T时刻，即最后计算的概率值
s11(T,1) = p1(T,1);
%由后往前计算，kim 的算法。
 for t=T-1:-1:1
     s00(t,1) = p0(t,1)*transP11*s00(t+1,1)/p00(t,1)...
         +p0(t,1)*(1-transP22)*s11(t+1,1)/p11(t,1);    
     s11(t,1) = p1(t,1)*(1-transP11)*s00(t+1,1)/p00(t,1)...
         +p1(t,1)*transP22*s11(t+1,1)/p11(t,1);    
 end

 smooth = [s00 s11];
 smoothed = smooth(2:end,:);
 %%
 %isreal(logL)

if isreal(logL) 
else
   logL = 1e6;
end

if isinf(logL)
   logL = 1e6;
end
if isnan(logL)
    logL = 1e6;
end