function [var_share] = vardcmp_wk(Q,sxi,D)
global T;global F;

nt=48;qnt=48*4;

Q1 = zeros(qnt,qnt);Q1(1,1)=Q(1,1);
Q2 = zeros(qnt,qnt);Q2(2,2)=Q(2,2);
Q3 = zeros(qnt,qnt);Q3(3,3)=Q(3,3);
Q4 = zeros(qnt,qnt);Q4(4,4)=Q(4,4);

for nn=1:nt
    tmp = 4*(nn-1);
    ssxi(tmp+1,nn:nn+T-1) = sxi(1,:);
    ssxi(tmp+2,nn:nn+T-1) = sxi(2,:);
    ssxi(tmp+3,nn:nn+T-1) = sxi(3,:);
    ssxi(tmp+4,nn:nn+T-1) = sxi(4,:);
end

 dssxi = ssxi(:,nt:T);
 dsxi = sxi(:,nt:T);
 TT = T-nt+1;
 nn = 49;
 Xi=dssxi;

 for i=1:nt
    xx(1,4*(i-1)+1) = i-1;
    xx(1,4*(i-1)+2) = i-1;
    xx(1,4*(i-1)+3) = i-1;
    xx(1,4*(i-1)+4) = i-1;
end

for i=1:qnt-1
tt(i+1,i:i+4*(nt)-1)=xx;
end

tau = tt(1:qnt,1:qnt);

vVY_total=zeros(nn,TT);
vVY_1=zeros(nn,TT);
vVY_2=zeros(nn,TT);
vVY_3=zeros(nn,TT);
vVY_2_q1=zeros(nn,TT);
vVY_2_q2=zeros(nn,TT);
vVY_2_q3=zeros(nn,TT);
vVY_2_q4=zeros(nn,TT);

for kk=1:48
Xq(:,:,kk)=(F^(kk-1))*Q*(F^(kk-1))';
Xq1(:,:,kk)=(F^(kk-1))*Q1*(F^(kk-1))';
Xq2(:,:,kk)=(F^(kk-1))*Q2*(F^(kk-1))';
Xq3(:,:,kk)=(F^(kk-1))*Q3*(F^(kk-1))';
Xq4(:,:,kk)=(F^(kk-1))*Q4*(F^(kk-1))';
end

for n=1:48
W1 = zeros(nt,TT);W2 = zeros(nt,TT);X1 = zeros(nt,TT);X2 = zeros(nt,TT);
MW1 = zeros(nt,nt,TT);MW2 = zeros(nt,nt,TT);
D_MW1 = zeros(nt,nt,TT);D_MW2 = zeros(nt,nt,TT);
MP1 = zeros(nt,nt,TT);MP2 = zeros(nt,nt,TT);
D_MP1 = zeros(nt,nt,TT);D_MP2 = zeros(nt,nt,TT);
MC1 = zeros(nt,nt,TT);MC2 = zeros(nt,nt,TT);
M_JH1 = zeros(nt,qnt,TT);
UL=zeros(nt,TT);UH=zeros(nt,TT);TU=zeros(nt,TT);
yhat = zeros(5,TT);
fJH1 = zeros(5,qnt,TT);fJH = zeros(qnt,5,TT);
SXq=zeros(qnt,qnt);SXq1=zeros(qnt,qnt);SXq2=zeros(qnt,qnt);SXq3=zeros(qnt,qnt);SXq4=zeros(qnt,qnt);

fVY_total=zeros(5,5,TT);
fVY_1=zeros(5,5,TT);
fVY_2=zeros(5,5,TT);
fVY_2_q1=zeros(5,5,TT);
fVY_2_q2=zeros(5,5,TT);
fVY_2_q3=zeros(5,5,TT);
fVY_2_q4=zeros(5,5,TT);
fVY_3=zeros(5,5,TT);

fxi0(:,n+1) = (F^n)*Xi(:,1);

for t=n+1:TT %TT

for i=1:nt
WL(i,t) = fxi0(1+4*(i-1),t);
WH(i,t) = fxi0(2+4*(i-1),t);
XL(i,t) = fxi0(3+4*(i-1),t);
XH(i,t) = fxi0(4+4*(i-1),t);
end

for i=1:nt
VXL(4*(i-1)+1,t)=XL(i,t);
VXL(4*(i-1)+2,t)=XL(i,t);
VXL(4*(i-1)+3,t)=XL(i,t);
VXL(4*(i-1)+4,t)=XL(i,t);
VXH(4*(i-1)+1,t)=XH(i,t);
VXH(4*(i-1)+2,t)=XH(i,t);
VXH(4*(i-1)+3,t)=XH(i,t);
VXH(4*(i-1)+4,t)=XH(i,t);
VWL(4*(i-1)+1,t)=WL(i,t);
VWL(4*(i-1)+2,t)=WL(i,t);
VWL(4*(i-1)+3,t)=WL(i,t);
VWL(4*(i-1)+4,t)=WL(i,t);
VWH(4*(i-1)+1,t)=WH(i,t);
VWH(4*(i-1)+2,t)=WH(i,t);
VWH(4*(i-1)+3,t)=WH(i,t);
VWH(4*(i-1)+4,t)=WH(i,t);
end

for i=1:qnt
    for j=1:qnt
        if i>j
      PML(i,j,t) = 1;
      PMH(i,j,t) = 1;
        else
        PML(i,j,t) = exp(-exp(VXL(i,t)+D(tau(i,j)+1,1)));
        PMH(i,j,t) = exp(-exp(VXH(i,t)+D(tau(i,j)+1,1)));
        end
    end
end

PML(1,:,t)=ones(1,qnt);
PMH(1,:,t)=ones(1,qnt);

% (4) computing weekly number unemployed: WU1, WU2 : in row "t", we have the weekly number unemployed at t%
% WU1, WU2: T by qnt matrix%

WUL(t,:) = VWL(:,t)'.*prod(PML(:,:,t),1); % weekly number unemployed of group 1
WUH(t,:) = VWH(:,t)'.*prod(PMH(:,:,t),1); % weekly number unemployed of group 2

%%%%%%%%%%%%
% jacobian %
%%%%%%%%%%%%

% jacobian wrt w1,w2 %

PPML(t,:) = prod(PML(:,:,t),1); % product of MPP1 matrix
PPMH(t,:) = prod(PMH(:,:,t),1); % product of MPP2 matrix

for i=1:nt
JPML(t,i) = PPML(t,4*(i-1)+1)+PPML(t,4*(i-1)+2)+PPML(t,4*(i-1)+3)+PPML(t,4*(i-1)+4);
JPMH(t,i) = PPMH(t,4*(i-1)+1)+PPMH(t,4*(i-1)+2)+PPMH(t,4*(i-1)+3)+PPMH(t,4*(i-1)+4);
end

% JH_W1/W2 matrix: rows indicate the index of measurement eqn, column indicate
% the variables %

JH_WL = zeros(nt,nt,T);
JH_WH = zeros(nt,nt,T);

for i=1:nt
JH_WL(i,i,t)=JPML(t,i);
JH_WH(i,i,t)=JPMH(t,i);
end

for i=1:qnt
    for j=1:qnt
        if i>j
      FML(i,j,t) = 0;
      FMH(i,j,t) = 0;
        else
        FML(i,j,t) = -exp(VXL(i,t)+D(tau(i,j)+1,1));
        FMH(i,j,t) = -exp(VXH(i,t)+D(tau(i,j)+1,1));
        end
    end
end

FML(1,:,t)=zeros(1,qnt);
FMH(1,:,t)=zeros(1,qnt);

% convert FML and FMH into nt by qnt matrices (sum 1~4 rows, 5~8 rows, ...
% etc) => MFML, MFMH

for i=1:nt
    MFML(i,:,t)=FML(4*(i-1)+1,:,t)+FML(4*(i-1)+2,:,t)+FML(4*(i-1)+3,:,t)+FML(4*(i-1)+4,:,t);
    MFMH(i,:,t)=FMH(4*(i-1)+1,:,t)+FMH(4*(i-1)+2,:,t)+FMH(4*(i-1)+3,:,t)+FMH(4*(i-1)+4,:,t);
end

% nt by qnt matrix
for i=1:nt
JW_XL(i,:,t) = MFML(i,:,t).*WUL(t,:);
JW_XH(i,:,t) = MFMH(i,:,t).*WUH(t,:);
end

% nt(variable) by nt(meas, number unemployed) matrix 

for j=1:nt
JM_XL(:,j,t) = JW_XL(:,4*(j-1)+1,t)+JW_XL(:,4*(j-1)+2,t)+JW_XL(:,4*(j-1)+3,t)+JW_XL(:,4*(j-1)+4,t); 
JM_XH(:,j,t) = JW_XH(:,4*(j-1)+1,t)+JW_XH(:,4*(j-1)+2,t)+JW_XH(:,4*(j-1)+3,t)+JW_XH(:,4*(j-1)+4,t); 
end

% now the row is monthly number unemployed, while column is variables %
JH_XL(:,:,t) = JM_XL(:,:,t)'; 
JH_XH(:,:,t) = JM_XH(:,:,t)';


for i=1:nt
    for j=1:nt
        M_JH1(i,4*(j-1)+1,t) = JH_WL(i,j,t);
        M_JH1(i,4*(j-1)+2,t) = JH_WH(i,j,t);
        M_JH1(i,4*(j-1)+3,t) = JH_XL(i,j,t);
        M_JH1(i,4*(j-1)+4,t) = JH_XH(i,j,t);
    end
end

fJH1(1,:,t)= M_JH1(1,:,t);
fJH1(2,:,t)= M_JH1(2,:,t)+M_JH1(3,:,t);
fJH1(3,:,t)= M_JH1(4,:,t)+M_JH1(5,:,t)+M_JH1(6,:,t);
fJH1(4,:,t)= sum(M_JH1(7:12,:,t),1);
fJH1(5,:,t)= sum(M_JH1(13:nt,:,t),1);

fJH(:,:,t)=fJH1(:,:,t)';

%(3)Variance of Y: E[(Y(t+m)-Y(t+m|t))(Y(t+m)-Y(t+m|t))'] 

SXq = sum(Xq(:,:,1:n),3);
SXq1 = sum(Xq1(:,:,1:n),3);
SXq2 = sum(Xq2(:,:,1:n),3);
SXq3 = sum(Xq3(:,:,1:n),3);
SXq4 = sum(Xq4(:,:,1:n),3);

fVY_2(:,:,t)=(fJH(:,:,t)')*(SXq)*(fJH(:,:,t));
fVY_2_q1(:,:,t)=(fJH(:,:,t)')*(SXq1)*(fJH(:,:,t));
fVY_2_q2(:,:,t)=(fJH(:,:,t)')*(SXq2)*(fJH(:,:,t));
fVY_2_q3(:,:,t)=(fJH(:,:,t)')*(SXq3)*(fJH(:,:,t));
fVY_2_q4(:,:,t)=(fJH(:,:,t)')*(SXq4)*(fJH(:,:,t));

vVY_2(n,t)=ones(1,5)*[(fJH(:,:,t)')*(SXq)*(fJH(:,:,t))]*ones(5,1);
vVY_2_q1(n,t)=ones(1,5)*[(fJH(:,:,t)')*(SXq1)*(fJH(:,:,t))]*ones(5,1);
vVY_2_q2(n,t)=ones(1,5)*[(fJH(:,:,t)')*(SXq2)*(fJH(:,:,t))]*ones(5,1);
vVY_2_q3(n,t)=ones(1,5)*[(fJH(:,:,t)')*(SXq3)*(fJH(:,:,t))]*ones(5,1);
vVY_2_q4(n,t)=ones(1,5)*[(fJH(:,:,t)')*(SXq4)*(fJH(:,:,t))]*ones(5,1);

fxi0(:,t+1) = (F^n)*Xi(:,t+1-n);
end

sum_vVY_q(n,1) = sum(vVY_2(n,:),2)/(TT-n);
sum_vVY_q(n,2) = sum(vVY_2_q1(n,:),2)/(TT-n);
sum_vVY_q(n,3) = sum(vVY_2_q2(n,:),2)/(TT-n);
sum_vVY_q(n,4) = sum(vVY_2_q3(n,:),2)/(TT-n);
sum_vVY_q(n,5) = sum(vVY_2_q4(n,:),2)/(TT-n);
end

for s=1:4
var_share(:,s) = sum_vVY_q(:,s+1)./sum_vVY_q(:,1);
end

end

