function [omegatilde,LL_RE,LL_ME,term1_RE_cum,LL_RE_cum,LL_ME_cum] = LL_vuong_MERE(datmat,N,T,xvec,xvec_RE,xvec_ME)

%Read in data
Y=datmat;
%Rearrange dimensions to get NxTx3-dimensional Y

%Minimize log likelihood 
rng(2501);
shocks.v=randn(1,1000);
shocks.w=randn(1,1000);

N=size(Y,1);    %number of forecasters
Nper=size(Y,2); %number of periods

%within forcaster-date settings:
Nz=3; 
Ns_RE=7; 
Ns_ME=6;

%reshape the data
fe=reshape(Y(:,:,1),[N*Nper,1]);
fr=reshape(Y(:,:,2),[N*Nper,1]);
lfe=reshape(Y(:,:,3),[N*Nper,1]);
yit=reshape([fe fr lfe]',[Nz*N,Nper]); 

%read in estimated parameters
mu=xvec(1);
rho1=xvec(2);
rho2=xvec(3);
sigw=xvec(4);
sigv=xvec_RE(1);
rhohat=xvec_ME(1);

%% Construct LL_RE
T=length(shocks.w);
s=zeros(1,T);
for t=3:T
    s(1,t) = rho1*s(1,t-1) + rho2*s(1,t-2) + sigw*shocks.w(1,t);
end
y = s + sigv*shocks.v;
AA = [rho1 rho2; 1 0]; 
BB = [sigw;0];
CC = [1 0];
DD = sigv;
s0=zeros(2,1); 
P0=eye(2);

% Run Kalman filter to get K_RE
t=1;
    %predict
    sbm1 = AA*s0;
    Pbm1 = AA*P0*AA' + BB*BB';
    %update
    z=y(1,t);
    K = Pbm1*CC'/(CC*Pbm1*CC'+DD*DD');
    sb = sbm1 + K*(z-CC*sbm1);
    Pb = (eye(2)-K*CC)*Pbm1;
for t=2:T
    %predict
     sbm1 = AA*sb;
     Pbm1 = AA*Pb*AA' + BB*BB';
       
    %update
    z=y(1,t);
    K = Pbm1*CC'/(CC*Pbm1*CC'+DD*DD');
    sb = sbm1 + K*(z-CC*sbm1);
    Pb = (eye(2)-K*CC)*Pbm1;
end 
k1_RE=K(1,1); 
k2_RE=K(2,1);


%% Run through filter to get K_ME
num = ((rho1-rhohat)^2+rho2^2)*(1-rho2) + 2*(rho1-rhohat)*rho2*rho1;
denom = (1-rho1^2-rho2^2)*(1-rho2) - 2*rho1^2*rho2;
Veps = sigw^2*(num/denom + 1);
sigeps=sqrt(Veps);

T=length(shocks.w);
s=zeros(1,T);
for t=3:T
    s(1,t) = rho1*s(1,t-1) + rho2*s(1,t-2) + sigw*shocks.w(1,t);
end
y=s + sigv*shocks.v;

AA = [rhohat 0; 1 0]; 
BB = [sigeps;0];
CC = [1 0];
DD = sigv;
s0=zeros(2,1); 
P0=eye(2);

t=1;
    %predict
    sbm1 = AA*s0;
    Pbm1 = AA*P0*AA' + BB*BB';
    %update
    z=y(1,t);
    K = Pbm1*CC'/(CC*Pbm1*CC'+DD*DD');
    sb = sbm1 + K*(z-CC*sbm1);
    Pb = (eye(2)-K*CC)*Pbm1;
for t=2:T
    %predict
     sbm1 = AA*sb;
     Pbm1 = AA*Pb*AA' + BB*BB';
       
    %update
    z=y(1,t);
    K = Pbm1*CC'/(CC*Pbm1*CC'+DD*DD');
    sb = sbm1 + K*(z-CC*sbm1);
    Pb = (eye(2)-K*CC)*Pbm1;
end 
k1_ME=K(1,1); 

%% Next run through another filter to construct LL_RE & LL_ME

% First set up matrices for RE
a_RE = [rho1 rho2 0 0 0 0 0;
    1 0 0 0 0 0 0;
    0 0 0 0 0 0 0;
    k1_RE*rho1 k1_RE*rho2 0 (1-k1_RE)*rho1 (1-k1_RE)*rho2 0 0;
    k2_RE*rho1 k2_RE*rho2 0 1-k2_RE*rho1 -k2_RE*rho2 0 0;
    0 0 0 rho1 rho2 0 0;
    0 0 0 1 0 0 0];

b_RE = [0 1 0;
    0 0 0;
    1 0 0;
    0 k1_RE k1_RE;
    0 k2_RE k2_RE;
    0 0 0;
    0 0 0];

c_RE = [rho1 rho2 1 -rho1 -rho2 0 0;
    0 0 0 rho1 rho2 -rho1 -rho2;
    0 1 0 0 0 0 -1];

%define relevant matrices (kron products are taken such that they are
%consistent with how yit is shaped)
T_RE = kron(speye(N),a_RE); 
B_RE = kron(speye(N),b_RE); 
Q_RE=kron(speye(N),[sigw^2 0 0; 0 sigw^2 0; 0 0 sigv^2]);
W_RE = (kron(speye(N),c_RE));

% Next set up matrices for ME
a_ME = [rho1 rho2 0 0 0 0;
    1 0 0 0 0 0;
    0 0 0 0 0 0;
    k1_ME*rho1 k1_ME*rho2 0 (1-k1_ME)*rhohat 0 0;
    0 0 0 rhohat 0 0;
    0 0 0 1 0 0];

b_ME = [0 1 0;
    0 0 0;
    1 0 0;
    0 k1_ME k1_ME;
    0 0 0;
    0 0 0];
c_ME = [rho1 rho2 1 -rhohat 0 0;
    0 0 0 rhohat -rhohat 0;
    0 1 0 0 0 -1];
%define relevant matrices (kron products are taken such that they are
%consistent with how yit is shaped)
T_ME = kron(speye(N),a_ME); 
B_ME = kron(speye(N),b_ME); 
Q_ME = kron(speye(N),[sigw^2 0 0; 0 sigw^2 0; 0 0 sigv^2]);
W_ME = (kron(speye(N),c_ME));

%% Now loop through observations and construct LL functions
omegatilde=0;
 
    %Setting initial state and variance 
    P0_RE=speye(N*Ns_RE); 
    F_RE = {}; 
    K_RE = {}; 
    Pttm1_RE = {}; 
    dttm1y_RE = {}; 
    d0y_RE = sparse(zeros(Ns_RE*N,1));  
    Syy_RE_cum = 0; 
    term1_RE_cum=0; 
    term2_RE_cum=0;

    P0_ME=speye(N*Ns_ME); 
    F_ME = {}; 
    K_ME = {}; 
    Pttm1_ME = {}; 
    dttm1y_ME = {}; 
    d0y_ME = sparse(zeros(Ns_ME*N,1));  
    Syy_ME_cum = 0; 
    term1_ME_cum=0; 
    term2_ME_cum=0;
    
   
    for t=1:Nper
        yt=yit(:,t);
        Mt=speye(Nz*N);
        Mt(isnan(yt),:)=[]; 
        yt(isnan(yt))=[];
        
        % RE filter
        Wt_RE=Mt*W_RE;
        if t==1
            F_RE{t} = Wt_RE*P0_RE*Wt_RE' ;
            K_RE{t} = T_RE*P0_RE*Wt_RE'/F_RE{t};
            ystar_RE{t} = yt - Wt_RE*d0y_RE;
            dttm1y_RE{t} = T_RE*d0y_RE + K_RE{t}*ystar_RE{t};
            Pttm1_RE{t} = T_RE*(P0_RE - P0_RE*Wt_RE'/F_RE{t}*Wt_RE*P0_RE)*T_RE' + B_RE*Q_RE*B_RE';
        else
            F_RE{t} = Wt_RE*Pttm1_RE{t-1}*Wt_RE' ;
            K_RE{t} = T_RE*Pttm1_RE{t-1}*Wt_RE'/F_RE{t};
            ystar_RE{t} = yt - Wt_RE*dttm1y_RE{t-1};
            dttm1y_RE{t} = T_RE*dttm1y_RE{t-1} + K_RE{t}*ystar_RE{t};
            Pttm1_RE{t} = T_RE*(Pttm1_RE{t-1} - Pttm1_RE{t-1}*Wt_RE'/F_RE{t}*Wt_RE*Pttm1_RE{t-1})*T_RE' + B_RE*Q_RE*B_RE';
        end
        
        % ME filter
        Wt_ME=Mt*W_ME;
        if t==1
            F_ME{t} = Wt_ME*P0_ME*Wt_ME' ;
            K_ME{t} = T_ME*P0_ME*Wt_ME'/F_ME{t};
            ystar_ME{t} = yt - Wt_ME*d0y_ME;
            dttm1y_ME{t} = T_ME*d0y_ME + K_ME{t}*ystar_ME{t};
            Pttm1_ME{t} = T_ME*(P0_ME - P0_ME*Wt_ME'/F_ME{t}*Wt_ME*P0_ME)*T_ME' + B_ME*Q_ME*B_ME';
        else
            F_ME{t} = Wt_ME*Pttm1_ME{t-1}*Wt_ME' ;
            K_ME{t} = T_ME*Pttm1_ME{t-1}*Wt_ME'/F_ME{t};
            ystar_ME{t} = yt - Wt_ME*dttm1y_ME{t-1};
            dttm1y_ME{t} = T_ME*dttm1y_ME{t-1} + K_ME{t}*ystar_ME{t};
            Pttm1_ME{t} = T_ME*(Pttm1_ME{t-1} - Pttm1_ME{t-1}*Wt_ME'/F_ME{t}*Wt_ME*Pttm1_ME{t-1})*T_ME' + B_ME*Q_ME*B_ME';
        end        
        
        % LL elements for RE
        term1_RE_tmp = size(Mt,1);
        term1_RE_cum = size(Mt,1) + term1_RE_cum;
        Syy_RE_tmp = ystar_RE{t}'/F_RE{t}*ystar_RE{t};
        Syy_RE_cum = ystar_RE{t}'/F_RE{t}*ystar_RE{t} + Syy_RE_cum;
        try
            term2_RE_tmp = 2*sum(log(diag(chol(F_RE{t})))) ;%determinant (Inf issues using det with large N)
            term2_RE_cum = 2*sum(log(diag(chol(F_RE{t})))) + term2_RE_cum;%determinant (Inf issues using det with large N)
        catch
            term2_RE_tmp = Inf;
            term2_RE_cum = Inf;
        end
        
        % LL elements for ME
        term1_ME_tmp = size(Mt,1);
        term1_ME_cum = size(Mt,1) + term1_ME_cum;
        Syy_ME_tmp = ystar_ME{t}'/F_ME{t}*ystar_ME{t};
        Syy_ME_cum = ystar_ME{t}'/F_ME{t}*ystar_ME{t} + Syy_ME_cum;
        try
            term2_ME_tmp = 2*sum(log(diag(chol(F_ME{t})))) ;%determinant (Inf issues using det with large N)
            term2_ME_cum = 2*sum(log(diag(chol(F_ME{t})))) + term2_ME_cum;%determinant (Inf issues using det with large N)
        catch
            term2_ME_tmp = Inf;
            term2_ME_cum = Inf;
        end
        
        % Compute LL_RE, LL_ME, and omegatilde
        LL_RE = -0.5 * ( term1_RE_tmp*log(2*pi) + term2_RE_tmp + Syy_RE_tmp) ;
        LL_ME = -0.5 * ( term1_ME_tmp*log(2*pi) + term2_ME_tmp + Syy_ME_tmp) ;
        omegatilde = (LL_ME-LL_RE)^2 + omegatilde;
    
    end %t
    
    omegatilde=1/term1_RE_cum*omegatilde;
    
    %Compute Log Likelihood
    LL_RE_cum = -0.5 * ( term1_RE_cum*log(2*pi) + term2_RE_cum + Syy_RE_cum) ;
    LL_ME_cum = -0.5 * ( term1_ME_cum*log(2*pi) + term2_ME_cum + Syy_ME_cum) ;
    

end
