function [omegatilde,LL_OC,LL_ME,term1_OC_cum,LL_OC_cum,LL_ME_cum] = LL_vuong_MEOC(datmat,N,T,xvec,xvec_RE,xvec_OC,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_OC=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);
alphav=xvec_OC(1);
rhohat=xvec_ME(1);

%% Construct LL_OC
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 = alphav*sigv;
s0=zeros(2,1); 
P0=eye(2);

% Run Kalman filter to get K_OC
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_OC=K(1,1); 
k2_OC=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_OC & LL_ME

% First set up matrices for OC
a_OC = [rho1 rho2 0 0 0 0 0;
    1 0 0 0 0 0 0;
    0 0 0 0 0 0 0;
    k1_OC*rho1 k1_OC*rho2 0 (1-k1_OC)*rho1 (1-k1_OC)*rho2 0 0;
    k2_OC*rho1 k2_OC*rho2 0 1-k2_OC*rho1 -k2_OC*rho2 0 0;
    0 0 0 rho1 rho2 0 0;
    0 0 0 1 0 0 0];

b_OC = [0 1 0;
    0 0 0;
    1 0 0;
    0 k1_OC k1_OC;
    0 k2_OC k2_OC;
    0 0 0;
    0 0 0];

c_OC = [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_OC = kron(speye(N),a_OC); 
B_OC = kron(speye(N),b_OC); 
Q_OC=kron(speye(N),[sigw^2 0 0; 0 sigw^2 0; 0 0 sigv^2]);
W_OC = (kron(speye(N),c_OC));

% 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_OC        = speye(N*Ns_OC); 
    F_OC         = {}; 
    K_OC         = {}; 
    Pttm1_OC     = {}; 
    dttm1y_OC    = {}; 
    d0y_OC       = sparse(zeros(Ns_OC*N,1));  
    Syy_OC_cum   = 0; 
    term1_OC_cum = 0; 
    term2_OC_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)) = [];
        
        % OC filter
        Wt_OC = Mt*W_OC;
        if t==1
            F_OC{t}      = Wt_OC*P0_OC*Wt_OC' ;
            K_OC{t}      = T_OC*P0_OC*Wt_OC'/F_OC{t};
            ystar_OC{t}  = yt - Wt_OC*d0y_OC;
            dttm1y_OC{t} = T_OC*d0y_OC + K_OC{t}*ystar_OC{t};
            Pttm1_OC{t}  = T_OC*(P0_OC - P0_OC*Wt_OC'/F_OC{t}*Wt_OC*P0_OC)*T_OC' + B_OC*Q_OC*B_OC';
        else
            F_OC{t}      = Wt_OC*Pttm1_OC{t-1}*Wt_OC' ;
            K_OC{t}      = T_OC*Pttm1_OC{t-1}*Wt_OC'/F_OC{t};
            ystar_OC{t}  = yt - Wt_OC*dttm1y_OC{t-1};
            dttm1y_OC{t} = T_OC*dttm1y_OC{t-1} + K_OC{t}*ystar_OC{t};
            Pttm1_OC{t}  = T_OC*(Pttm1_OC{t-1} - Pttm1_OC{t-1}*Wt_OC'/F_OC{t}*Wt_OC*Pttm1_OC{t-1})*T_OC' + B_OC*Q_OC*B_OC';
        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 OC
        term1_OC_tmp = size(Mt,1);
        term1_OC_cum = size(Mt,1) + term1_OC_cum;
        Syy_OC_tmp   = ystar_OC{t}'/F_OC{t}*ystar_OC{t};
        Syy_OC_cum   = ystar_OC{t}'/F_OC{t}*ystar_OC{t} + Syy_OC_cum;
        try
            term2_OC_tmp = 2*sum(log(diag(chol(F_OC{t})))) ;%determinant (Inf issues using det with large N)
            term2_OC_cum = 2*sum(log(diag(chol(F_OC{t})))) + term2_OC_cum;%determinant (Inf issues using det with large N)
        catch
            term2_OC_tmp = Inf;
            term2_OC_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_OC, LL_ME, and omegatilde
        LL_OC = -0.5 * ( term1_OC_tmp*log(2*pi) + term2_OC_tmp + Syy_OC_tmp) ;
        LL_ME = -0.5 * ( term1_ME_tmp*log(2*pi) + term2_ME_tmp + Syy_ME_tmp) ;
        omegatilde = (LL_ME-LL_OC)^2 + omegatilde;
    
    end %t
    
    omegatilde=1/term1_OC_cum*omegatilde;
    
    %Compute Log Likelihood
    LL_OC_cum = -0.5 * ( term1_OC_cum*log(2*pi) + term2_OC_cum + Syy_OC_cum) ;
    LL_ME_cum = -0.5 * ( term1_ME_cum*log(2*pi) + term2_ME_cum + Syy_ME_cum) ;
    

end
