%% 1. ESTIMATORS: MAIN FUNCTIONS

%% 1.1 Second-order spatial lag simultaneous equations bivariate probit 
function [b1S,b2S,g1S,g2S,l11S,l12S,l21S,l22S,tauS,ma1S,ma2S,vcva11S,vcva12S,vcva22S,y1starS,y2starS,alpha1S,alpha2S] = sho_sem_biprobit_2W(y1,y2,X1,X2,W1,W2,id,year,iter) 

    % Second-order spatial lag simultaneous equations bivariate probit 

    % INPUT
    % y1,y2: Binary dependent variable for equation 1 and 2
    % X1,X2: Regressors for equation 1 and 2
    % W1,W2: Spatial weight matrices 
    % id: individual observation
    % year: year
    % iter: number of MCMC iterations
    
    % OUTPUT
    % b1S,b2S: MCMC draws for beta1 and beta2
    % g1S,g2S: MCMC draws for gamma1 and gamma2
    % l11S,l12S: MCMC draws for lambda11 and lambda12 (equation 1)
    % l21S,l22S: MCMC draws for lambda21 and lambda22 (equation 2)
    % tauS: MCMC draws for tau
    % alpha1S,alphpa2S: MCMC draws for alpha1 and alpha2
    % y1starS,y2starS: MCMC draws for y1* and y2*
    % ma1S,ma2S: averages of MCMC draws for alpha1 and alpha2
    % vcva11S,vcva12S,vcva22S: covariances of MCMC draws for alpha1
    % and alpha2
  
    nsample=1;

    % Data preparations
    un=size(unique(id),1);                    
    nt=size(year,1);                          
    year=year-min(year)+1;
    t=max(year);
    n_per_t = arrayfun(@(i) sum(year==i),(1:t)','UniformOutput',false); 
    Ic=cellfun(@(Y) eye(Y),n_per_t,'UniformOutput',false);
    npyr_end=cumsum(cell2mat(n_per_t));
    npyr_start=npyr_end-cell2mat(n_per_t)+1;
    temp=[npyr_start npyr_end]; 
    npy_range = arrayfun(@(i) temp,(1:t)','UniformOutput',false);
    npy_range = npy_range{1,1};

    clear npyr_end npyr_start
    
    y1c = arrayfun(@(i) y1(npy_range(i,1):npy_range(i,2),:),...
        (1:t)','UniformOutput',false);   
    y2c = arrayfun(@(i) y2(npy_range(i,1):npy_range(i,2),:),...
        (1:t)','UniformOutput',false);
    X1c = arrayfun(@(i) X1(npy_range(i,1):npy_range(i,2),:),...
        (1:t)','UniformOutput',false);
    X2c = arrayfun(@(i) X2(npy_range(i,1):npy_range(i,2),:),...
        (1:t)','UniformOutput',false);
    W1c = arrayfun(@(i) W1(npy_range(i,1):npy_range(i,2),npy_range(i,1):npy_range(i,2)),...
        (1:t)','UniformOutput',false);  
    W2c = arrayfun(@(i) W2(npy_range(i,1):npy_range(i,2),npy_range(i,1):npy_range(i,2)),...
        (1:t)','UniformOutput',false);

    y1ind=y1c;
    y2ind=y2c;

    A = zeros(nt,un);
    uniqueid = unique(id);

    for i = 1:un
        I = find(id == uniqueid(i));
        A(I,i) = 1;
    end              

    yperun=sum(A,1);

    % Preallocate baskets for results
    [tauS,g1S,g2S,ma1S,ma2S,vcva11S,vcva12S,vcva22S,...
        l11S,l12S,l21S,l22S]=deal(zeros(iter,1,'single'));

    [y1starS,y2starS,alpha1S,alpha2S]=deal(zeros(nt,iter,'single'));

    k1=size(X1,2); 
    k2=size(X2,2); 
    b1S=zeros(iter,k1,'single');
    b2S=zeros(iter,k2,'single');


    % Priors
    c1 = zeros(k1,1); V_b10=eye(k1)*1e+12; V_b10I=inv(V_b10); V_b10Ic1=V_b10I*c1;
    c2 = zeros(k2,1); V_b20=eye(k2)*1e+12; V_b20I=inv(V_b20); V_b20Ic2=V_b20I*c2;
    mu_mua10=0; V_mua10=1;                    
    mu_mua20=0; V_mua20=1;
    mu_mua0=[mu_mua10;mu_mua20];
    V_mua0=[V_mua10 0; 0 V_mua20];
    V_mua0I=inv(V_mua0);
    V_alphaI0=eye(2);
    v_Va0=2;   

    % Start values and vectors
    l11=0.1; ccl11=0.1; accl11=0;
    l12=0.1; ccl12=0.1; accl12=0;
    l21=0.1; ccl21=0.1; accl21=0;
    l22=0.1; ccl22=0.1; accl22=0; 
    tau=0.2; cctau=0.5; acctau=0;
    g1=0.1; ccg1=0.1; accg1=0;             
    g2=-0.1; ccg2=0.1; accg2=0;
    
    alpha=zeros(2*un,1); 
    V_alphaI=eye(2);

    [z1,z2]=deal(zeros(nt,1)); 
    z1c=arrayfun(@(i) z1(npy_range(i,1):npy_range(i,2),:),...
            (1:t)','UniformOutput',false);
    z2c=arrayfun(@(i) z2(npy_range(i,1):npy_range(i,2),:),...
            (1:t)','UniformOutput',false);       

    clear z1 z2

    % Dependent matrices and values
    psi11=(1-tau^2)^(-1);
    psi22=psi11;
    psi12=-tau*psi11;
    psi21=psi12;
    psitil=[psi11 psi12; psi21 psi22];

    alpha_nt1=A*alpha(1:un,1);
    alpha_nt1c = arrayfun(@(i) alpha_nt1(npy_range(i,1):npy_range(i,2),:),...
        (1:t)','UniformOutput',false);
    alpha_nt2=A*alpha(un+1:2*un,1);
    alpha_nt2c = arrayfun(@(i) alpha_nt2(npy_range(i,1):npy_range(i,2),:),...
        (1:t)','UniformOutput',false);

    Vb=zeros(k1+k2,k1+k2);
    F=zeros(2,2);

    X1TX1=X1'*X1; X1TX2=X1'*X2; X2TX1=X2'*X1; X2TX2=X2'*X2;
    VA=eye(2*un,2*un);

    L11 = cellfun(@(A,B,C) eye(A)-l11*B-l12*C,n_per_t,W1c,W2c,...
        'UniformOutput',false);
    L22 = cellfun(@(A,B,C) eye(A)-l21*B-l22*C,n_per_t,W1c,W2c,...
        'UniformOutput',false);
    
    W1W1 = cellfun(@(W) W*W,W1c,'UniformOutput',false);
    W2W2 = cellfun(@(W) W*W,W2c,'UniformOutput',false);
    W1W2 = cellfun(@(U,W) U*W,W1c,W2c,'UniformOutput',false);
    W2W1 = cellfun(@(U,W) U*W,W2c,W1c,'UniformOutput',false);

    LYAc1=cellfun(@(U,V,W,Z) U*V-g1*W-Z, L11,y1c,y2c, alpha_nt1c,...
        'UniformOutput', false);
    LYAc2=cellfun(@(U,V,W,Z)-g2*U+V*W-Z, y1c,L22,y2c, alpha_nt2c, ...
        'UniformOutput', false);  

    clear  indus* uniqueid y1 y2

    nW1=norm(W1);
    nW2=norm(W2);

    % MCMC ESTIMATION
    for m=1:iter

        % Parameters on regressors
        Vb(1:k1,1:k1)=psi11*X1TX1+V_b10I;
        Vb(1:k1,k1+1:end)=psi12*X1TX2;
        Vb(k1+1:end,1:k1)=psi21*X2TX1;
        Vb(k1+1:end,k1+1:end)=psi22*X2TX2;+V_b10I;

        Vb=(Vb+Vb')/2;
        Vb=inv(Vb);

        mutemp1=cellfun(@(U,W,Y) U'*(psi11*W+psi12*Y), X1c,LYAc1,LYAc2, 'UniformOutput', false); 
        mutemp2=cellfun(@(U,W,Y) U'*(psi21*W+psi22*Y), X2c,LYAc1,LYAc2, 'UniformOutput', false); 
        matSize = size(mutemp1{1},1);
        B = reshape(cell2mat(mutemp1'),matSize,[],t);
        mubet1=sum(B,3)+V_b10Ic1;
        matSize = size(mutemp2{1},1);
        B = reshape(cell2mat(mutemp2'),matSize,[],t);
        mubet2=sum(B,3)+V_b20Ic2;    
        mubeta=Vb*([mubet1; mubet2]);

        beta=mvnrnd(mubeta,Vb);

        beta1=beta(1,1:k1)';
        beta2=beta(1,k1+1:end)';
        b1S(m,:)=beta1';   
        b2S(m,:)=beta2';
        
        L11L22=cellfun(@(U,V,W,X,Z,Y,A) U-(l11+l21)*V-(l12+l22)*W+l11*l21*X...
            +l12*l21*Z+l22*l11*Y+l12*l22*A,Ic,W1c,W2c,W1W1,W2W1,W1W2,W2W2,'UniformOutput', false); 

        f1=cellfun(@(U,V) U-V*beta1, LYAc1,X1c, 'UniformOutput', false);
        f2=cellfun(@(U,V) U-V*beta2, LYAc2,X2c, 'UniformOutput', false);
        f1f1=cellfun(@(U) U'*U, f1, 'UniformOutput', false);
        f1f2=cellfun(@(U,V) U'*V, f1,f2, 'UniformOutput', false);
        f2f2=cellfun(@(U) U'*U, f2, 'UniformOutput', false);
        F(1,1)=sum(cell2mat(f1f1));
        F(1,2)=sum(cell2mat(f1f2));
        F(2,1)=F(1,2);
        F(2,2)=sum(cell2mat(f2f2));

        expterm=0.5*(trace(F*psitil));

        tempD=cellfun(@(U,V) U-g1*g2*V, L11L22,Ic, 'UniformOutput', false); 
        detm =  arrayfun(@(i) det(tempD{i,1}),(1:t)','UniformOutput',false);
        detm = prod(cell2mat(detm));
        detm = log(detm);  
        cd1 =  detm - expterm;

        % Gammas
        cond1=(1+abs(l11)*nW1+abs(l12)*nW2)*(1+abs(l21)*nW1+abs(l22)*nW2);
        cond2=(1-abs(l11)*nW1-abs(l12)*nW2)*(1-abs(l21)*nW1-abs(l22)*nW2);
        
        ftemp1=cellfun(@(U,V) U+g1*V, f1,y2c, 'UniformOutput', false);
        [g1,ccg1,accg1,expterm,f1,F] = s_p_g1(g1,g2,cond1,cond2,y2c,f1,f2,ftemp1,psitil,ccg1,accg1,expterm,F,m,t,Ic,L11L22,cd1);
        g1S(m,1)=g1;

        ftemp2=cellfun(@(U,V) U+g2*V, f2,y1c, 'UniformOutput', false);  
        [g2,ccg2,accg2,f2,expterm,F,cd1] = s_p_g2(g2,g1,cond1,cond2,y1c,f1,f2,ftemp2,psitil,ccg2,accg2,expterm,F,m,t,Ic,L11L22,cd1);
        g2S(m,1)=g2;
        
        % Lambdas
        g1g2Ic=cellfun(@(V) g1*g2*V, Ic,'UniformOutput', false);
        W1y1=cellfun(@(U,V) U*V, W1c,y1c,'UniformOutput', false);
        W2y1=cellfun(@(U,V) U*V, W2c,y1c,'UniformOutput', false);
        W1y2=cellfun(@(U,V) U*V, W1c,y2c,'UniformOutput', false);
        W2y2=cellfun(@(U,V) U*V, W2c,y2c,'UniformOutput', false);

        L11y1=cellfun(@(V,W) V*W, L11,y1c, 'UniformOutput', false); 

        pr=1-abs(l12);
        ftemp1=cellfun(@(U,V) U-V, f1,L11y1, 'UniformOutput', false);
        [l11,L11,L11y1,L11L22,ccl11,accl11,expterm,f1,F,cd1] = s_p_l11_2W(l11,l12,pr,L11y1,L11,L11L22,y1c,f1,f2,ftemp1,psitil,W1c,W2c,ccl11,accl11,expterm,F,m,t,cd1,Ic,g1g2Ic,W1W1,W2W1,W1W2,W2W2,l21,l22,W1y1,W2y1);
        l11S(m,1)=l11; 

        pr=1-abs(l11);
        ftemp1=cellfun(@(U,V) U-V, f1,L11y1, 'UniformOutput', false);
        [l12,L11,L11y1,L11L22,ccl12,accl12,expterm,f1,F,cd1] = s_p_l12_2W(l12,l11,pr,L11y1,L11,L11L22,y1c,f1,f2,ftemp1,psitil,W1c,W2c,ccl12,accl12,expterm,F,m,t,cd1,Ic,g1g2Ic,W1W1,W2W1,W1W2,W2W2,l21,l22,W1y1,W2y1);
        l12S(m,1)=l12;
        
        L22y2=cellfun(@(V,W) V*W, L22,y2c, 'UniformOutput', false); 

        pr=1-abs(l22);
        ftemp2=cellfun(@(U,V) U-V, f2,L22y2, 'UniformOutput', false);
        [l21,L22,L22y2,L11L22,ccl21,accl21,expterm,f2,F,cd1] = s_p_l21_2W(l21,l22,pr,L22y2,L22,L11L22,y2c,f1,f2,ftemp2,psitil,W1c,W2c,ccl21,accl21,expterm,F,m,t,cd1,Ic,g1g2Ic,W1W1,W2W1,W1W2,W2W2,l11,l12,W1y2,W2y2);
        l21S(m,1)=l21; 

        pr=1-abs(l21);
        ftemp2=cellfun(@(U,V) U-V, f2,L22y2, 'UniformOutput', false);
        [l22,L22,L22y2,~,ccl22,accl22,expterm,~,F,~] = s_p_l22_2W(l22,l21,pr,L22y2,L22,L11L22,y2c,f1,f2,ftemp2,psitil,W1c,W2c,ccl22,accl22,expterm,F,m,t,cd1,Ic,g1g2Ic,W1W1,W2W1,W1W2,W2W2,l11,l12,W1y2,W2y2);
        l22S(m,1)=l22; 

        % Tau
        pr1=1;
        [tau,cctau,acctau,psitil]=s_p_tau(tau,pr1,F,psitil,cctau,acctau,expterm,nt,m);
        tauS(m,1)=tau; 
        psi11=psitil(1,1);
        psi12=psitil(1,2);
        psi21=psitil(2,1);
        psi22=psitil(2,2);

        YXB1=cellfun(@(V,W,Z)  V-g1*W-Z*beta1, L11y1, y2c, X1c, 'UniformOutput', false);
        YXB2=cellfun(@(V,W,Z) -g2*V+W-Z*beta2, y1c, L22y2, X2c, 'UniformOutput', false);

        % Alphas
        mtemp1=cellfun(@(V,Z) psi11*V+psi12*Z, YXB1,YXB2, 'UniformOutput', false);  
        mtemp2=cellfun(@(V,Z) psi21*V+psi22*Z, YXB1,YXB2, 'UniformOutput', false);
        mtemp1=A'*cell2mat(mtemp1);
        mtemp2=A'*cell2mat(mtemp2);

        [alpha,V_alphaI] = h_alpha(alpha(1:un,1),alpha(un+1:2*un,1),...
            psitil,mtemp1,mtemp2,V_alphaI,mu_mua0,V_mua0I,v_Va0,V_alphaI0,un,yperun,VA);

        ma1S(m,:)=mean(alpha(1:un,1));
        ma2S(m,:)=mean(alpha(un+1:end,1));
        vcvtemp=cov(alpha(1:un,1),alpha(un+1:end,1));
        vcva11S(m,:)=vcvtemp(1,1);
        vcva22S(m,:)=vcvtemp(2,2);
        vcva12S(m,:)=vcvtemp(1,2);

        alpha_nt1=A*alpha(1:un,1);
        alpha_nt2=A*alpha(un+1:end,1);

        alpha_nt1c = arrayfun(@(i) alpha_nt1(npy_range(i,1):npy_range(i,2),:),(1:t)','UniformOutput',false);
        alpha_nt2c = arrayfun(@(i) alpha_nt2(npy_range(i,1):npy_range(i,2),:),(1:t)','UniformOutput',false);
        alpha1S(:,m)=cell2mat(alpha_nt1c);
        alpha2S(:,m)=cell2mat(alpha_nt2c);

        % Ystars
        XBA1=cellfun(@(U,V) U*beta1+V,X1c,alpha_nt1c,'UniformOutput', false);
        XBA2=cellfun(@(U,V) U*beta2+V,X2c,alpha_nt2c,'UniformOutput', false);

        L11I=cellfun(@(V)  inv(V), L11, 'UniformOutput', false);
        Ltil22=cellfun(@(U,V) inv(U-g1*g2*V), L22, L11I, 'UniformOutput', false);
        Ltil21=cellfun(@(U,V) g2*U*V, Ltil22, L11I, 'UniformOutput', false);
        Ltil12=cellfun(@(U,V) g1*U*V, L11I, Ltil22, 'UniformOutput', false);
        Ltil11=cellfun(@(U,V) U+g1*U*V,L11I,Ltil21,'UniformOutput', false);

        Ltil_1111=cellfun(@(U) U*U', Ltil11,'UniformOutput', false);  
        Ltil_1211=cellfun(@(U,V) U*V', Ltil12,Ltil11,'UniformOutput', false);
        Ltil_1212=cellfun(@(U) U*U', Ltil12,'UniformOutput', false);  
        V11=cellfun(@(U,V,W) U+tau*(V+V')+W, Ltil_1111,Ltil_1211,Ltil_1212,'UniformOutput', false);  

        Ltil_1121=cellfun(@(U,V) U*V', Ltil11, Ltil21,'UniformOutput', false);  
        Ltil_1221=cellfun(@(U,V) U*V', Ltil12, Ltil21,'UniformOutput', false);  
        Ltil_1122=cellfun(@(U,V) U*V', Ltil11, Ltil22,'UniformOutput', false);
        Ltil_1222=cellfun(@(U,V) U*V', Ltil12, Ltil22,'UniformOutput', false);
        V12=cellfun(@(U,V,W,Z) U+tau*(V+W)+Z, Ltil_1121,Ltil_1221,Ltil_1122,Ltil_1222,'UniformOutput', false);  
        V21=cellfun(@(U) U', V12,'UniformOutput', false); 

        Ltil_2121=cellfun(@(U) U*U', Ltil21,'UniformOutput', false);
        Ltil_2221=cellfun(@(U,V) U*V', Ltil22, Ltil21,'UniformOutput', false);
        Ltil_2222=cellfun(@(U) U*U', Ltil22,'UniformOutput', false);
        V22=cellfun(@(U,V,W) U+tau*(V+V')+W, Ltil_2121,Ltil_2221,Ltil_2222,'UniformOutput', false); 

        mu1=cellfun(@(U,V,W,Z) U*V+W*Z, Ltil11,XBA1,Ltil12,XBA2,'UniformOutput', false); 
        mu2=cellfun(@(U,V,W,Z) U*V+W*Z, Ltil21,XBA1,Ltil22,XBA2,'UniformOutput', false); 

        muy1=cellfun(@(U,V,W,X,Y) U+V*(W\(X-Y)),mu1,V12,V22,y2c,mu2,'UniformOutput', false);
        Vy1=cellfun(@(U,V,W,X) U-V*inv(W)*X, V11,V12,V22,V21,'UniformOutput', false); 
        Vy1=cellfun(@(U) inv(U), Vy1,'UniformOutput', false); 
        [y1c,z1c]=arrayfun(@(i) draw_ystar(y1ind{i,1},z1c{i,1},muy1{i,1},Vy1{i,1},n_per_t{i,1},nsample),(1:t)','UniformOutput',false);
        y1starS(:,m)=cell2mat(y1c);

        muy2=cellfun(@(U,V,W,X,Y) U+V*(W\(X-Y)),mu2,V21,V11,y1c,mu1,'UniformOutput', false);
        Vy2=cellfun(@(U,V,W,X) U-V*inv(W)*X, V22,V21,V11,V12,'UniformOutput', false); 
        Vy2=cellfun(@(U) inv(U), Vy2,'UniformOutput', false); 
        [y2c,z2c]=arrayfun(@(i) draw_ystar(y2ind{i,1},z2c{i,1},muy2{i,1},Vy2{i,1},n_per_t{i,1},nsample),(1:t)','UniformOutput',false);
        y2starS(:,m)=cell2mat(y2c);

        LYAc1=cellfun(@(U,V,W,Z) U*V-g1*W-Z, L11,y1c,y2c, alpha_nt1c, 'UniformOutput', false);
        LYAc2=cellfun(@(U,V,W,Z)-g2*U+V*W-Z, y1c,L22,y2c, alpha_nt2c, 'UniformOutput', false);  

    end

end



%% 1.2 Third-order spatial lag simultaneous equations bivariate probit 
function [b1S,b2S,g1S,g2S,l11S,l12S,l13S,l21S,l22S,l23S,tauS,ma1S,ma2S,vcva11S,vcva12S,vcva22S,y1starS,y2starS,alpha1S,alpha2S] = sho_sem_biprobit_3W(y1,y2,X1,X2,W1,W2,W3,id,year,iter) 

    % Third-order spatial lag simultaneous equations bivariate probit 
    
    % INPUT
    % y1,y2: Binary dependent variable for equation 1 and 2
    % X1,X2: Regressors for equation 1 and 2
    % W1,W2,W3: Spatial weight matrices 
    % id: individual observation
    % year: year
    % iter: number of MCMC iterations
    
    % OUTPUT
    % b1S,b2S: MCMC draws for beta1 and beta2
    % g1S,g2S: MCMC draws for gamma1 and gamma2
    % l11S,l12S,l13S: MCMC draws for lambda11, lambda12, and lambda13 (equation 1)
    % l21S,l22S,l23S: MCMC draws for lambda21, lambda22, and lambda23 (equation 2)
    % tauS: MCMC draws for tau
    % alpha1S,alphpa2S: MCMC draws for alpha1 and alpha2
    % y1starS,y2starS: MCMC draws for y1* and y2*
    % ma1S,ma2S: averages of MCMC draws for alpha1 and alpha2
    % vcva11S,vcva12S,vcva22S: covariances of MCMC draws for alpha1
    % and alpha2
  
    nsample=1;

    % Data preparations
    un=size(unique(id),1);                    
    nt=size(year,1);                          
    year=year-min(year)+1;
    t=max(year);

    n_per_t = arrayfun(@(i) sum(year==i),(1:t)','UniformOutput',false); 

    Ic=cellfun(@(Y) eye(Y),n_per_t,'UniformOutput',false);

    npyr_end=cumsum(cell2mat(n_per_t));
    npyr_start=npyr_end-cell2mat(n_per_t)+1;
    temp=[npyr_start npyr_end]; 
    npy_range = arrayfun(@(i) temp,(1:t)','UniformOutput',false);
    npy_range = npy_range{1,1};

    clear npyr_end npyr_start

    y1c = arrayfun(@(i) y1(npy_range(i,1):npy_range(i,2),:),...
        (1:t)','UniformOutput',false);   
    y2c = arrayfun(@(i) y2(npy_range(i,1):npy_range(i,2),:),...
        (1:t)','UniformOutput',false);
    X1c = arrayfun(@(i) X1(npy_range(i,1):npy_range(i,2),:),...
        (1:t)','UniformOutput',false);
    X2c = arrayfun(@(i) X2(npy_range(i,1):npy_range(i,2),:),...
        (1:t)','UniformOutput',false);
    W1c = arrayfun(@(i) W1(npy_range(i,1):npy_range(i,2),npy_range(i,1):npy_range(i,2)),...
        (1:t)','UniformOutput',false);  
    W2c = arrayfun(@(i) W2(npy_range(i,1):npy_range(i,2),npy_range(i,1):npy_range(i,2)),...
        (1:t)','UniformOutput',false);
    W3c = arrayfun(@(i) W3(npy_range(i,1):npy_range(i,2),npy_range(i,1):npy_range(i,2)),...
        (1:t)','UniformOutput',false);
    
    y1ind=y1c;
    y2ind=y2c;

    A = zeros(nt,un);
    uniqueid = unique(id);

    for i = 1:un
        I = find(id == uniqueid(i));
        A(I,i) = 1;
    end              

    yperun=sum(A,1);

    % Preallocate the baskets for results
    [tauS,g1S,g2S,ma1S,ma2S,vcva11S,vcva12S,vcva22S,...
        l11S,l12S,l13S,l21S,l22S,l23S]=deal(zeros(iter,1,'single'));

    [y1starS,y2starS,alpha1S,alpha2S]=deal(zeros(nt,iter,'single'));

    k1=size(X1,2); 
    k2=size(X2,2); 
    b1S=zeros(iter,k1,'single');
    b2S=zeros(iter,k2,'single');

    % Priors
    c1 = zeros(k1,1); V_b10=eye(k1)*1e+12; V_b10I=inv(V_b10); V_b10Ic1=V_b10I*c1;
    c2 = zeros(k2,1); V_b20=eye(k2)*1e+12; V_b20I=inv(V_b20); V_b20Ic2=V_b20I*c2;
    mu_mua10=0; V_mua10=1;                    
    mu_mua20=0; V_mua20=1;
    mu_mua0=[mu_mua10;mu_mua20];
    V_mua0=[V_mua10 0; 0 V_mua20];
    V_mua0I=inv(V_mua0);
    V_alphaI0=eye(2);
    v_Va0=2;   

    % Start values and vectors
    l11=0.1; ccl11=0.1; accl11=0;
    l12=0.1; ccl12=0.1; accl12=0;
    l13=0.1; ccl13=0.1; accl13=0;
    l21=0.1; ccl21=0.1; accl21=0;
    l22=0.1; ccl22=0.1; accl22=0; 
    l23=0.1; ccl23=0.1; accl23=0; 
    tau=0.2; cctau=0.5; acctau=0;
    g1=0.1; ccg1=0.1; accg1=0;             
    g2=-0.1; ccg2=0.1; accg2=0;
    
    alpha=zeros(2*un,1); 
    V_alphaI=eye(2);

    [z1,z2]=deal(zeros(nt,1)); 
    z1c=arrayfun(@(i) z1(npy_range(i,1):npy_range(i,2),:),...
            (1:t)','UniformOutput',false);
    z2c=arrayfun(@(i) z2(npy_range(i,1):npy_range(i,2),:),...
            (1:t)','UniformOutput',false);       

    clear z1 z2

    % Dependent matrices and values
    psi11=(1-tau^2)^(-1);
    psi22=psi11;
    psi12=-tau*psi11;
    psi21=psi12;
    psitil=[psi11 psi12; psi21 psi22];

    alpha_nt1=A*alpha(1:un,1);
    alpha_nt1c = arrayfun(@(i) alpha_nt1(npy_range(i,1):npy_range(i,2),:),...
        (1:t)','UniformOutput',false);
    alpha_nt2=A*alpha(un+1:2*un,1);
    alpha_nt2c = arrayfun(@(i) alpha_nt2(npy_range(i,1):npy_range(i,2),:),...
        (1:t)','UniformOutput',false);

    Vb=zeros(k1+k2,k1+k2);
    F=zeros(2,2);

    X1TX1=X1'*X1; X1TX2=X1'*X2; X2TX1=X2'*X1; X2TX2=X2'*X2;
    VA=eye(2*un,2*un);

    L11 = cellfun(@(A,B,C,D) eye(A)-l11*B-l12*C-l13*D,n_per_t,W1c,W2c,W3c,...
        'UniformOutput',false);
    L22 = cellfun(@(A,B,C,D) eye(A)-l21*B-l22*C-l23*D,n_per_t,W1c,W2c,W3c,...
        'UniformOutput',false);
    
    W1W1 = cellfun(@(W) W*W,W1c,'UniformOutput',false);
    W2W2 = cellfun(@(W) W*W,W2c,'UniformOutput',false);
    W3W3 = cellfun(@(W) W*W,W3c,'UniformOutput',false);
    W1W2 = cellfun(@(U,W) U*W,W1c,W2c,'UniformOutput',false);
    W2W1 = cellfun(@(U,W) U*W,W2c,W1c,'UniformOutput',false);
    W1W3 = cellfun(@(U,W) U*W,W1c,W3c,'UniformOutput',false);
    W3W1 = cellfun(@(U,W) U*W,W3c,W1c,'UniformOutput',false);
    W2W3 = cellfun(@(U,W) U*W,W2c,W3c,'UniformOutput',false);
    W3W2 = cellfun(@(U,W) U*W,W3c,W2c,'UniformOutput',false);

    LYAc1=cellfun(@(U,V,W,Z) U*V-g1*W-Z, L11,y1c,y2c, alpha_nt1c,...
        'UniformOutput', false);
    LYAc2=cellfun(@(U,V,W,Z)-g2*U+V*W-Z, y1c,L22,y2c, alpha_nt2c, ...
        'UniformOutput', false);  

    clear  indus* uniqueid y1 y2

    nW1=norm(W1);
    nW2=norm(W2);
    nW3=norm(W3);

    % MCMC ESTIMATION
    for m=1:iter

        % Parameters on regressors
        Vb(1:k1,1:k1)=psi11*X1TX1+V_b10I;
        Vb(1:k1,k1+1:end)=psi12*X1TX2;
        Vb(k1+1:end,1:k1)=psi21*X2TX1;
        Vb(k1+1:end,k1+1:end)=psi22*X2TX2;+V_b10I;
        Vb=(Vb+Vb')/2;
        Vb=inv(Vb);
        mutemp1=cellfun(@(U,W,Y) U'*(psi11*W+psi12*Y), X1c,LYAc1,LYAc2, 'UniformOutput', false); 
        mutemp2=cellfun(@(U,W,Y) U'*(psi21*W+psi22*Y), X2c,LYAc1,LYAc2, 'UniformOutput', false); 
        matSize = size(mutemp1{1},1);
        B = reshape(cell2mat(mutemp1'),matSize,[],t);
        mubet1=sum(B,3)+V_b10Ic1;
        matSize = size(mutemp2{1},1);
        B = reshape(cell2mat(mutemp2'),matSize,[],t);
        mubet2=sum(B,3)+V_b20Ic2;    
        mubeta=Vb*([mubet1; mubet2]);

        beta=mvnrnd(mubeta,Vb);

        beta1=beta(1,1:k1)';
        beta2=beta(1,k1+1:end)';
        b1S(m,:)=beta1';   
        b2S(m,:)=beta2';

        L11L22=cellfun(@(A,B,C,D,E,F,G,H,I,K,L,M,N) ...
            A-(l11+l21)*B-(l12+l22)*C-(l13+l23)*D ...
            +l11*l21*E+l11*l22*F+l11*l23*G...
            +l12*l21*H+l12*l22*I+l12*l23*K...
            +l13*l21*L+l13*l22*M+l13*l23*N,...
            Ic,W1c,W2c,W3c,W1W1,W1W2,W1W3,W2W1,W2W2,W2W3,W3W1,W3W2,W3W3,...
            'UniformOutput', false); 

        f1=cellfun(@(U,V) U-V*beta1, LYAc1,X1c, 'UniformOutput', false);
        f2=cellfun(@(U,V) U-V*beta2, LYAc2,X2c, 'UniformOutput', false);
        f1f1=cellfun(@(U) U'*U, f1, 'UniformOutput', false);
        f1f2=cellfun(@(U,V) U'*V, f1,f2, 'UniformOutput', false);
        f2f2=cellfun(@(U) U'*U, f2, 'UniformOutput', false);
        F(1,1)=sum(cell2mat(f1f1));
        F(1,2)=sum(cell2mat(f1f2));
        F(2,1)=F(1,2);
        F(2,2)=sum(cell2mat(f2f2));

        expterm=0.5*(trace(F*psitil));

        tempD=cellfun(@(U,V) U-g1*g2*V, L11L22,Ic, 'UniformOutput', false); 
        detm =  arrayfun(@(i) det(tempD{i,1}),(1:t)','UniformOutput',false);
        detm = prod(cell2mat(detm));
        detm = log(detm);  
        cd1 =  detm - expterm;

        % Gammas
        cond1=(1+abs(l11)*nW1+abs(l12)*nW2+abs(l13)*nW3)*(1+abs(l21)*nW1+abs(l22)*nW2+abs(l23)*nW3);
        cond2=(1-abs(l11)*nW1-abs(l12)*nW2-abs(l13)*nW3)*(1-abs(l21)*nW1-abs(l22)*nW2-abs(l23)*nW3);

        ftemp1=cellfun(@(U,V) U+g1*V, f1,y2c, 'UniformOutput', false);
        [g1,ccg1,accg1,expterm,f1,F] = s_p_g1(g1,g2,cond1,cond2,y2c,f1,f2,ftemp1,psitil,ccg1,accg1,expterm,F,m,t,Ic,L11L22,cd1);
        g1S(m,1)=g1;

        ftemp2=cellfun(@(U,V) U+g2*V, f2,y1c, 'UniformOutput', false);  
        [g2,ccg2,accg2,f2,expterm,F,cd1] = s_p_g2(g2,g1,cond1,cond2,y1c,f1,f2,ftemp2,psitil,ccg2,accg2,expterm,F,m,t,Ic,L11L22,cd1);
        g2S(m,1)=g2; 

        % Lambdas
        g1g2Ic=cellfun(@(V) g1*g2*V, Ic,'UniformOutput', false);
        W1y1=cellfun(@(U,V) U*V, W1c,y1c,'UniformOutput', false);
        W2y1=cellfun(@(U,V) U*V, W2c,y1c,'UniformOutput', false);
        W3y1=cellfun(@(U,V) U*V, W3c,y1c,'UniformOutput', false);
        W1y2=cellfun(@(U,V) U*V, W1c,y2c,'UniformOutput', false);
        W2y2=cellfun(@(U,V) U*V, W2c,y2c,'UniformOutput', false);
        W3y2=cellfun(@(U,V) U*V, W3c,y2c,'UniformOutput', false);
        L11y1=cellfun(@(V,W) V*W, L11,y1c, 'UniformOutput', false); 

        pr=1-abs(l12)-abs(l13);
        ftemp1=cellfun(@(U,V) U-V, f1,L11y1, 'UniformOutput', false);
        [l11,L11,L11y1,L11L22,ccl11,accl11,expterm,f1,F,cd1] = s_p_l11_3W(l11,l12,l13,pr,L11y1,L11,L11L22,y1c,f1,f2,ftemp1,psitil,ccl11,accl11,expterm,F,m,t,cd1,Ic,g1g2Ic,W1c,W2c,W3c,W1W1,W1W2,W1W3,W2W1,W2W2,W2W3,W3W1,W3W2,W3W3,l21,l22,l23,W1y1,W2y1,W3y1);
        l11S(m,1)=l11; 

        pr=1-abs(l11)-abs(l13);
        ftemp1=cellfun(@(U,V) U-V, f1,L11y1, 'UniformOutput', false);
        [l12,L11,L11y1,L11L22,ccl12,accl12,expterm,f1,F,cd1] = s_p_l12_3W(l12,l11,l13,pr,L11y1,L11,L11L22,y1c,f1,f2,ftemp1,psitil,ccl12,accl12,expterm,F,m,t,cd1,Ic,g1g2Ic,W1c,W2c,W3c,W1W1,W1W2,W1W3,W2W1,W2W2,W2W3,W3W1,W3W2,W3W3,l21,l22,l23,W1y1,W2y1,W3y1);
        l12S(m,1)=l12;

        pr=1-abs(l11)-abs(l12);
        ftemp1=cellfun(@(U,V) U-V, f1,L11y1, 'UniformOutput', false);
        [l13,L11,L11y1,L11L22,ccl13,accl13,expterm,f1,F,cd1] = s_p_l13_3W(l13,l11,l12,pr,L11y1,L11,L11L22,y1c,f1,f2,ftemp1,psitil,ccl13,accl13,expterm,F,m,t,cd1,Ic,g1g2Ic,W1c,W2c,W3c,W1W1,W1W2,W1W3,W2W1,W2W2,W2W3,W3W1,W3W2,W3W3,l21,l22,l23,W1y1,W2y1,W3y1);
        l13S(m,1)=l13; 

        L22y2=cellfun(@(V,W) V*W, L22,y2c, 'UniformOutput', false); 

        pr=1-abs(l22)-abs(l23);
        ftemp2=cellfun(@(U,V) U-V, f2,L22y2, 'UniformOutput', false);
        [l21,L22,L22y2,L11L22,ccl21,accl21,expterm,f2,F,cd1] = s_p_l21_3W(l21,l22,l23,pr,L22y2,L22,L11L22,y2c,f1,f2,ftemp2,psitil,ccl21,accl21,expterm,F,m,t,cd1,Ic,g1g2Ic,W1c,W2c,W3c,W1W1,W1W2,W1W3,W2W1,W2W2,W2W3,W3W1,W3W2,W3W3,l11,l12,l13,W1y2,W2y2,W3y2);
        l21S(m,1)=l21; 

        pr=1-abs(l21)-abs(l23);
        ftemp2=cellfun(@(U,V) U-V, f2,L22y2, 'UniformOutput', false);
        [l22,L22,L22y2,L11L22,ccl22,accl22,expterm,f2,F,cd1] = s_p_l22_3W(l22,l21,l23,pr,L22y2,L22,L11L22,y2c,f1,f2,ftemp2,psitil,ccl22,accl22,expterm,F,m,t,cd1,Ic,g1g2Ic,W1c,W2c,W3c,W1W1,W1W2,W1W3,W2W1,W2W2,W2W3,W3W1,W3W2,W3W3,l11,l12,l13,W1y2,W2y2,W3y2);
        l22S(m,1)=l22; 

        pr=1-abs(l21)-abs(l22);
        ftemp2=cellfun(@(U,V) U-V, f2,L22y2, 'UniformOutput', false);
        [l23,L22,L22y2,~,ccl23,accl23,expterm,~,F,~] = s_p_l23_3W(l23,l21,l22,pr,L22y2,L22,L11L22,y2c,f1,f2,ftemp2,psitil,ccl23,accl23,expterm,F,m,t,cd1,Ic,g1g2Ic,W1c,W2c,W3c,W1W1,W1W2,W1W3,W2W1,W2W2,W2W3,W3W1,W3W2,W3W3,l11,l12,l13,W1y2,W2y2,W3y2);
        l23S(m,1)=l23; 

        % Tau
        pr1=1;
        [tau,cctau,acctau,psitil]=s_p_tau(tau,pr1,F,psitil,cctau,acctau,expterm,nt,m);

        tauS(m,1)=tau; 
        psi11=psitil(1,1);
        psi12=psitil(1,2);
        psi21=psitil(2,1);
        psi22=psitil(2,2);

        YXB1=cellfun(@(V,W,Z)  V-g1*W-Z*beta1, L11y1, y2c, X1c, 'UniformOutput', false);
        YXB2=cellfun(@(V,W,Z) -g2*V+W-Z*beta2, y1c, L22y2, X2c, 'UniformOutput', false);

        % Alphas
        mtemp1=cellfun(@(V,Z) psi11*V+psi12*Z, YXB1,YXB2, 'UniformOutput', false);  
        mtemp2=cellfun(@(V,Z) psi21*V+psi22*Z, YXB1,YXB2, 'UniformOutput', false);
        mtemp1=A'*cell2mat(mtemp1);
        mtemp2=A'*cell2mat(mtemp2);

        [alpha,V_alphaI] = h_alpha(alpha(1:un,1),alpha(un+1:2*un,1),...
            psitil,mtemp1,mtemp2,V_alphaI,mu_mua0,V_mua0I,v_Va0,V_alphaI0,un,yperun,VA);

        ma1S(m,:)=mean(alpha(1:un,1));
        ma2S(m,:)=mean(alpha(un+1:end,1));
        vcvtemp=cov(alpha(1:un,1),alpha(un+1:end,1));
        vcva11S(m,:)=vcvtemp(1,1);
        vcva22S(m,:)=vcvtemp(2,2);
        vcva12S(m,:)=vcvtemp(1,2);

        alpha_nt1=A*alpha(1:un,1);
        alpha_nt2=A*alpha(un+1:end,1);

        alpha_nt1c = arrayfun(@(i) alpha_nt1(npy_range(i,1):npy_range(i,2),:),(1:t)','UniformOutput',false);
        alpha_nt2c = arrayfun(@(i) alpha_nt2(npy_range(i,1):npy_range(i,2),:),(1:t)','UniformOutput',false);
        alpha1S(:,m)=cell2mat(alpha_nt1c);
        alpha2S(:,m)=cell2mat(alpha_nt2c);
        
        % Ystars
        XBA1=cellfun(@(U,V) U*beta1+V,X1c,alpha_nt1c,'UniformOutput', false);
        XBA2=cellfun(@(U,V) U*beta2+V,X2c,alpha_nt2c,'UniformOutput', false);
        L11I=cellfun(@(V)  inv(V), L11, 'UniformOutput', false);
        Ltil22=cellfun(@(U,V) inv(U-g1*g2*V), L22, L11I, 'UniformOutput', false);
        Ltil21=cellfun(@(U,V) g2*U*V, Ltil22, L11I, 'UniformOutput', false);
        Ltil12=cellfun(@(U,V) g1*U*V, L11I, Ltil22, 'UniformOutput', false);
        Ltil11=cellfun(@(U,V) U+g1*U*V,L11I,Ltil21,'UniformOutput', false);

        Ltil_1111=cellfun(@(U) U*U', Ltil11,'UniformOutput', false);  
        Ltil_1211=cellfun(@(U,V) U*V', Ltil12,Ltil11,'UniformOutput', false);
        Ltil_1212=cellfun(@(U) U*U', Ltil12,'UniformOutput', false);  
        V11=cellfun(@(U,V,W) U+tau*(V+V')+W, Ltil_1111,Ltil_1211,Ltil_1212,'UniformOutput', false);  

        Ltil_1121=cellfun(@(U,V) U*V', Ltil11, Ltil21,'UniformOutput', false);  
        Ltil_1221=cellfun(@(U,V) U*V', Ltil12, Ltil21,'UniformOutput', false);  
        Ltil_1122=cellfun(@(U,V) U*V', Ltil11, Ltil22,'UniformOutput', false);
        Ltil_1222=cellfun(@(U,V) U*V', Ltil12, Ltil22,'UniformOutput', false);
        V12=cellfun(@(U,V,W,Z) U+tau*(V+W)+Z, Ltil_1121,Ltil_1221,Ltil_1122,Ltil_1222,'UniformOutput', false);  
        V21=cellfun(@(U) U', V12,'UniformOutput', false); 

        Ltil_2121=cellfun(@(U) U*U', Ltil21,'UniformOutput', false);
        Ltil_2221=cellfun(@(U,V) U*V', Ltil22, Ltil21,'UniformOutput', false);
        Ltil_2222=cellfun(@(U) U*U', Ltil22,'UniformOutput', false);
        V22=cellfun(@(U,V,W) U+tau*(V+V')+W, Ltil_2121,Ltil_2221,Ltil_2222,'UniformOutput', false); 

        mu1=cellfun(@(U,V,W,Z) U*V+W*Z, Ltil11,XBA1,Ltil12,XBA2,'UniformOutput', false); 
        mu2=cellfun(@(U,V,W,Z) U*V+W*Z, Ltil21,XBA1,Ltil22,XBA2,'UniformOutput', false); 

        muy1=cellfun(@(U,V,W,X,Y) U+V*(W\(X-Y)),mu1,V12,V22,y2c,mu2,'UniformOutput', false);
        Vy1=cellfun(@(U,V,W,X) U-V*inv(W)*X, V11,V12,V22,V21,'UniformOutput', false); 
        Vy1=cellfun(@(U) inv(U), Vy1,'UniformOutput', false); 
        [y1c,z1c]=arrayfun(@(i) draw_ystar(y1ind{i,1},z1c{i,1},muy1{i,1},Vy1{i,1},n_per_t{i,1},nsample),(1:t)','UniformOutput',false);
        y1starS(:,m)=cell2mat(y1c);

        muy2=cellfun(@(U,V,W,X,Y) U+V*(W\(X-Y)),mu2,V21,V11,y1c,mu1,'UniformOutput', false);
        Vy2=cellfun(@(U,V,W,X) U-V*inv(W)*X, V22,V21,V11,V12,'UniformOutput', false); 
        Vy2=cellfun(@(U) inv(U), Vy2,'UniformOutput', false); 
        [y2c,z2c]=arrayfun(@(i) draw_ystar(y2ind{i,1},z2c{i,1},muy2{i,1},Vy2{i,1},n_per_t{i,1},nsample),(1:t)','UniformOutput',false);
        y2starS(:,m)=cell2mat(y2c);

        LYAc1=cellfun(@(U,V,W,Z) U*V-g1*W-Z, L11,y1c,y2c, alpha_nt1c, 'UniformOutput', false);
        LYAc2=cellfun(@(U,V,W,Z)-g2*U+V*W-Z, y1c,L22,y2c, alpha_nt2c, 'UniformOutput', false);  

    end

end



%% 2. SUBFUNCTIONS

%% 2.1 Same subfunctions for both models

function [g1,ccl,accl,expterm,f1,F,cd1] = s_p_g1(g1,g2,cond1,cond2,y2,f1,f2,ftemp1,psitil,ccl,accl,expterm,F,m,t,Ic,L11L22,cd1)
    
    % Metropolis Hastings procedure for gamma1

    g1c = g1+ ccl*randn(1,1);
    
    if (abs(g1c)*abs(g2) < cond1) && (abs(g1c)*abs(g2) > cond2)
        p=0;   
    else    
        f1c=cellfun(@(U,V) -g1c*U+V, y2,ftemp1, 'UniformOutput', false); 
        f1f1=cellfun(@(U) U'*U, f1c, 'UniformOutput', false);
        f1f2=cellfun(@(U,V) U'*V, f1c,f2, 'UniformOutput', false);
        Fc=F;
        Fc(1,1)=sum(cell2mat(f1f1));
        Fc(1,2)=sum(cell2mat(f1f2));
        Fc(2,1)=Fc(1,2);
        exptermc=0.5*(trace(Fc*psitil));
        tempD=cellfun(@(U,V) U-g1c*g2*V, L11L22,Ic, 'UniformOutput', false); 
        detm =  arrayfun(@(i) det(tempD{i,1}),(1:t)','UniformOutput',false);
        detm = log(prod(cell2mat(detm)));  
        cd1c =  detm - exptermc;
        p = exp(cd1c-cd1);
        p = min(1,p);

    end
    
    if p>unifrnd(0,1)
        g1 = g1c;
        cd1=cd1c;
        accl = accl+1; 
        f1=f1c;
        expterm=exptermc;
        F=Fc;
    end

    
    if accl/m < 0.4
          ccl = ccl/1.1;
    elseif accl/m > 0.6
          ccl = ccl*1.1;
    end   
   
end


function [g2,ccl,accl,f2,expterm,F,cd1] = s_p_g2(g2,g1,cond1,cond2,y1,f1,f2,ftemp2,psitil,ccl,accl,expterm,F,m,t,Ic,L11L22,cd1)

    % Metropolis Hastings procedure for gamma2
    
    g2c = g2+ ccl*randn(1,1);
 
    if (abs(g1)*abs(g2c) < cond1) && (abs(g1)*abs(g2c) > cond2)
        p=0;  
    else    
        f2c=cellfun(@(U,V) -g2c*U+V, y1,ftemp2, 'UniformOutput', false);  
        f2f2=cellfun(@(U) U'*U, f2c, 'UniformOutput', false);
        f1f2=cellfun(@(U,V) U'*V, f1,f2c, 'UniformOutput', false);
        Fc=F;
        Fc(2,2)=sum(cell2mat(f2f2));
        Fc(1,2)=sum(cell2mat(f1f2));
        Fc(2,1)=Fc(1,2);
        exptermc=0.5*(trace(Fc*psitil));
        tempD=cellfun(@(U,V) U-g1*g2c*V, L11L22,Ic, 'UniformOutput', false); 
        detm =  arrayfun(@(i) det(tempD{i,1}),(1:t)','UniformOutput',false);
        detm = log(prod(cell2mat(detm)));  
        cd1c =  detm - exptermc;
        p = exp(cd1c-cd1);
        p = min(1,p);
    end
    
    if  unifrnd(0,1)<p
        g2 = g2c;
        cd1 = cd1c;
        accl = accl+1; 
        f2=f2c;
        expterm=exptermc;
        F=Fc;   
    end
    
    if accl/m < 0.4
          ccl = ccl/1.1;
    elseif accl/m > 0.6
          ccl = ccl*1.1;
    end    
   
end



function [tau,ccl,accl,psitil] = s_p_tau(tau,pr1,F,psitil,ccl,accl,expterm,ntot,m)
    
    %  Metropolis Hastings procedure for tau
    
    tauc = tau+ ccl*randn(1,1);

    if abs(tauc) > pr1
        p=0;      
    else
        detm = -(ntot/2)*log(1-tau^2);
        cd1 =  detm - expterm;
        psic=[1 tauc; tauc 1];
        psitilc=inv(psic);
        exptermc=0.5*(trace(F*psitilc));
        detmc = -(ntot/2)*log(1-tauc^2);
        cd1c =  detmc - exptermc;
        p = exp(cd1c-cd1);
        p = min(1,p);  
    end
   
    if p>unifrnd(0,1)
        tau = tauc;
        accl = accl+1; 
        psitil=psitilc;
   
    end

    if accl/m < 0.4
          ccl = ccl/1.1;
    elseif accl/m > 0.6
          ccl = ccl*1.1;
    elseif accl/m<0.2
          ccl = 0.5;
    end   

end



function [alpha,V_alphaI] = h_alpha(alpha1,alpha2,psitil,mutemp1,mutemp2,V_alphaI,mu_mua0,V_mua0I,v0,V_alphaI0,un,yperid,VA)

    % Alpha draw (hierarchical structure)

    twun=2*un;
    
    % mu_alpha
    V_mua=inv(un*V_alphaI+V_mua0I);
    
    temp=[V_alphaI(1,1)*sum(alpha1)+V_alphaI(1,2)*sum(alpha2); ...
          V_alphaI(2,1)*sum(alpha1)+V_alphaI(2,2)*sum(alpha2)];
    mu_mua=V_mua*(temp+V_mua0I*mu_mua0);
    
    mu_al=mvnrnd(mu_mua,V_mua);
    
    % V_alpha
    v1=v0+un;
    b1=bsxfun(@minus,alpha1,mu_al(1,1));
    b2=bsxfun(@minus,alpha2,mu_al(1,2));
   
    H = [b1'*b1 b1'*b2;...
         b2'*b1 b2'*b2 ];
   
    temp=H+V_alphaI0;
    temp=(temp+temp')/2;
    VVal=inv(temp);

    V_alphaI=wishrnd(VVal,v1);
   
    % Alpha
    VA(1:un,1:un)=diag(psitil(1,1)*yperid+V_alphaI(1,1));
    VA(1:un,un+1:end)=diag(psitil(1,2)*yperid+V_alphaI(1,2));
    VA(un+1:end,1:un)=diag(psitil(2,1)*yperid+V_alphaI(2,1));
    VA(un+1:end,un+1:end)=diag(psitil(2,2)*yperid+V_alphaI(2,2));

    mutemp1=mutemp1+(V_alphaI(1,:)*mu_al')*ones(un,1);
    mutemp2=mutemp2+(V_alphaI(2,:)*mu_al')*ones(un,1);
    
    R=chol(VA);
    mua=R\(R'\[mutemp1;mutemp2]);
   
    tempa=randn(twun,1);
    alpha=R\tempa+mua;
   
end



function [y,z] = draw_ystar(yin,z,mu,Psi,n,nsample)

    % Draw ystar
    aa=diag(Psi);
    h=1./sqrt(aa);
    aatemp = - 1./aa;
    c=bsxfun(@times,Psi,aatemp);
    c(1:n+1:n*n) = 0;

    for nsiter=1:nsample
        
        gaz=c*z;
        zo1=z;
        
        for i = 1:n
            muuse = (-mu(i,1)-gaz(i,1))/h(i,1);
            m_left=-6*(1-yin(i,1))+yin(i,1).*muuse;
            m_right=(1-yin(i,1)).*muuse+6*yin(i,1);  
            a= 0.5*(1+erf(m_left/sqrt(2)));       
            b= 0.5*(1+erf(m_right/sqrt(2)));      
            p = a + (b-a)*rand(1);
            v = sqrt(2) * erfinv(2*p-1);           

            if isinf(v)==0 && isnan(v)==0
                z(i,1) = gaz(i,1) + h(i,1)*v; 
                gaz=gaz+c(:,i)*(z(i)-zo1(i));
            end      
                 
        end
       
    end
    
    y=mu+z;
    
end


%% 2.2 Lambda draws for second-order spatial lag simultaneous equations bivariate probit

function [l11,L11,L11y1,L11L22,ccl,accl,expterm,f1,F,cd1] = s_p_l11_2W(l11,l12,pr1,L11y1,L11,L11L22,y1,f1,f2,ftemp1,psitil,W1,W2,ccl,accl,expterm,F,m,t,cd1,Ic,g1g2Ic,W1W1,W2W1,W1W2,W2W2,l21,l22,W1y1,W2y1)
     
    % Metropolis Hastings procedure for lambda 11 
    
    l11c = l11+ ccl*randn(1,1);
  
    if abs(l11c) > pr1     
        p=0;          
    else 
        L11y1c=cellfun(@(U,V,W) U-l11c*V-l12*W, y1,W1y1,W2y1, 'UniformOutput', false);
        f1c=cellfun(@(U,V) U+V, L11y1c,ftemp1, 'UniformOutput', false); 
        f1f1=cellfun(@(U) U'*U, f1c, 'UniformOutput', false);
        f1f2=cellfun(@(U,V) U'*V, f1c,f2, 'UniformOutput', false);
        Fc=F;
        Fc(1,1)=sum(cell2mat(f1f1));
        Fc(1,2)=sum(cell2mat(f1f2));
        Fc(2,1)=Fc(1,2);
        exptermc=0.5*(trace(Fc*psitil));
        L11L22c=cellfun(@(U,V,W,X,Z,Y,A) U-(l11c+l21)*V-(l12+l22)*W+l11c*l21*X...
            +l12*l21*Z+l22*l11c*Y+l12*l22*A,Ic,W1,W2,W1W1,W2W1,W1W2,W2W2,'UniformOutput', false); 
        tempD=cellfun(@(U,V) U-V, L11L22c,g1g2Ic, 'UniformOutput', false); 
        detm =  arrayfun(@(i) det(tempD{i,1}),(1:t)','UniformOutput',false);
        detm = log(prod(cell2mat(detm)));  
        cd1c =  detm - exptermc;
        p = exp(cd1c-cd1);
        p = min(1,p);
    end
   

    if p>unifrnd(0,1)
        dl=l11-l11c; 
        l11 = l11c; 
        L11=cellfun(@(U,V) U+dl*V, L11,W1,'UniformOutput',false);
        L11y1=L11y1c;
        accl = accl+1; 
        expterm=exptermc;
        f1=f1c;
        F=Fc;
        cd1=cd1c;
        L11L22=L11L22c;
    end
    
    if accl/m < 0.4
          ccl = ccl/1.1;
    elseif accl/m > 0.6
          ccl = ccl*1.1;
    end  

end


function [l12,L11,L11y1,L11L22,ccl,accl,expterm,f1,F,cd1] = s_p_l12_2W(l12,l11,pr1,L11y1,L11,L11L22,y1,f1,f2,ftemp1,psitil,W1,W2,ccl,accl,expterm,F,m,t,cd1,Ic,g1g2Ic,W1W1,W2W1,W1W2,W2W2,l21,l22,W1y1,W2y1)
     
    % Metropolis Hastings procedure for lambda 12 
    l12c = l12+ ccl*randn(1,1);
  
    if abs(l12c) > pr1
        p=0;      
    else 
        L11y1c=cellfun(@(U,V,W) U-l11*V-l12c*W, y1,W1y1,W2y1, 'UniformOutput', false);
        f1c=cellfun(@(U,V) U+V, L11y1c,ftemp1, 'UniformOutput', false); 
        f1f1=cellfun(@(U) U'*U, f1c, 'UniformOutput', false);
        f1f2=cellfun(@(U,V) U'*V, f1c,f2, 'UniformOutput', false);
        Fc=F;
        Fc(1,1)=sum(cell2mat(f1f1));
        Fc(1,2)=sum(cell2mat(f1f2));
        Fc(2,1)=Fc(1,2);
        exptermc=0.5*(trace(Fc*psitil));
        L11L22c=cellfun(@(U,V,W,X,Z,Y,A) U-(l11+l21)*V-(l12c+l22)*W+l11*l21*X...
            +l12c*l21*Z+l22*l11*Y+l12c*l22*A,Ic,W1,W2,W1W1,W2W1,W1W2,W2W2,'UniformOutput', false); 
        tempD=cellfun(@(U,V) U-V, L11L22c,g1g2Ic, 'UniformOutput', false); 
        detm =  arrayfun(@(i) det(tempD{i,1}),(1:t)','UniformOutput',false);
        detm = log(prod(cell2mat(detm)));  
        cd1c =  detm - exptermc;
        p = exp(cd1c-cd1);
        p = min(1,p);
    end
   
    
    if p>unifrnd(0,1)
        dl=l12-l12c; 
        L11=cellfun(@(U,V) U+dl*V, L11,W2,'UniformOutput',false);
        l12 = l12c;
        L11y1=L11y1c;
        accl = accl+1; 
        expterm=exptermc;
        f1=f1c;
        F=Fc;
        cd1=cd1c;
        L11L22=L11L22c;
    end
    
    if accl/m < 0.4
          ccl = ccl/1.1;
    elseif accl/m > 0.6
          ccl = ccl*1.1;
    end  

end


function [l21,L22,L22y2,L11L22,ccl,accl,expterm,f2,F,cd1] = s_p_l21_2W(l21,l22,pr1,L22y2,L22,L11L22,y2,f1,f2,ftemp2,psitil,W1,W2,ccl,accl,expterm,F,m,t,cd1,Ic,g1g2Ic,W1W1,W2W1,W1W2,W2W2,l11,l12,W1y2,W2y2)
    
    % Metropolis Hastings procedure for lambda 21 
    l21c = l21+ ccl*randn(1,1);
    
    if abs(l21c) > pr1
        p=0;     
    else
        L22y2c=cellfun(@(U,V,W) U-l21c*V-l22*W, y2,W1y2,W2y2, 'UniformOutput', false);
        f2c=cellfun(@(U,V) U+V, L22y2c,ftemp2, 'UniformOutput', false); 
        f1f2=cellfun(@(U,V) U'*V, f1,f2c, 'UniformOutput', false);
        f2f2=cellfun(@(U) U'*U, f2c, 'UniformOutput', false); 
        Fc=F;
        Fc(1,2)=sum(cell2mat(f1f2));
        Fc(2,1)=Fc(1,2);
        Fc(2,2)=sum(cell2mat(f2f2));
        exptermc=0.5*(trace(Fc*psitil));
        L11L22c=cellfun(@(U,V,W,X,Z,Y,A) U-(l11+l21c)*V-(l12+l22)*W+l11*l21c*X...
            +l12*l21c*Z+l22*l11*Y+l12*l22*A,Ic,W1,W2,W1W1,W2W1,W1W2,W2W2,'UniformOutput', false); 
        tempD=cellfun(@(U,V) U-V, L11L22c,g1g2Ic, 'UniformOutput', false); 
        detm =  arrayfun(@(i) det(tempD{i,1}),(1:t)','UniformOutput',false);
        detm =  log(prod(cell2mat(detm)));  
        cd1c =  detm - exptermc;
        p = exp(cd1c-cd1);
        p = min(1,p);
    end
   
    if p>unifrnd(0,1)
        dl=l21-l21c; 
        l21 = l21c; 
        L22=cellfun(@(U,V) U+dl*V, L22,W1,'UniformOutput',false);
        L22y2=L22y2c;
        L11L22=L11L22c;
        accl = accl+1; 
        expterm=exptermc;
        f2=f2c;
        F=Fc;
        cd1=cd1c;
    end

    if accl/m < 0.4
          ccl = ccl/1.1;
    elseif accl/m > 0.6
          ccl = ccl*1.1;
    end  

end


function [l22,L22,L22y2,L11L22,ccl,accl,expterm,f2,F,cd1] = s_p_l22_2W(l22,l21,pr1,L22y2,L22,L11L22,y2,f1,f2,ftemp2,psitil,W1,W2,ccl,accl,expterm,F,m,t,cd1,Ic,g1g2Ic,W1W1,W2W1,W1W2,W2W2,l11,l12,W1y2,W2y2)
    
    % Metropolis Hastings procedure for lambda 22 
    l22c = l22+ ccl*randn(1,1);
    
    if abs(l22c) > pr1
        p=0;     
    else
        L22y2c=cellfun(@(U,V,W) U-l21*V-l22c*W, y2,W1y2,W2y2, 'UniformOutput', false);
        f2c=cellfun(@(U,V) U+V, L22y2c,ftemp2, 'UniformOutput', false); 
        f1f2=cellfun(@(U,V) U'*V, f1,f2c, 'UniformOutput', false);
        f2f2=cellfun(@(U) U'*U, f2c, 'UniformOutput', false); 
        Fc=F;
        Fc(1,2)=sum(cell2mat(f1f2));
        Fc(2,1)=Fc(1,2);
        Fc(2,2)=sum(cell2mat(f2f2));
        exptermc=0.5*(trace(Fc*psitil));
        L11L22c=cellfun(@(U,V,W,X,Z,Y,A) U-(l11+l21)*V-(l12+l22c)*W+l11*l21*X...
            +l12*l21*Z+l22c*l11*Y+l12*l22c*A,Ic,W1,W2,W1W1,W2W1,W1W2,W2W2,'UniformOutput', false); 
        tempD=cellfun(@(U,V) U-V, L11L22c,g1g2Ic, 'UniformOutput', false); 
        detm =  arrayfun(@(i) det(tempD{i,1}),(1:t)','UniformOutput',false);
        detm =  log(prod(cell2mat(detm)));  
        cd1c =  detm - exptermc;
        p = exp(cd1c-cd1);
        p = min(1,p);
    end
   
    if p>unifrnd(0,1)
        dl=l22-l22c; 
        l22 = l22c; 
        L22=cellfun(@(U,V) U+dl*V, L22,W2,'UniformOutput',false);
        L22y2=L22y2c;
        L11L22=L11L22c;
        accl = accl+1; 
        expterm=exptermc;
        f2=f2c;
        F=Fc;
        cd1=cd1c;
    end

    if accl/m < 0.4
          ccl = ccl/1.1;
    elseif accl/m > 0.6
          ccl = ccl*1.1;
    end  

end



%% 2.3 Lambda draws for third-order spatial lag simultaneous equations bivariate probit

function [l11,L11,L11y1,L11L22,ccl,accl,expterm,f1,F,cd1] = s_p_l11_3W(l11,l12,l13,pr1,L11y1,L11,L11L22,y1,f1,f2,ftemp1,psitil,ccl,accl,expterm,F,m,t,cd1,Ic,g1g2Ic,W1,W2,W3,W1W1,W1W2,W1W3,W2W1,W2W2,W2W3,W3W1,W3W2,W3W3,l21,l22,l23,W1y1,W2y1,W3y1)
     
    % Metropolis Hastings procedure for lambda 11
    
    l11c = l11+ ccl*randn(1,1);
  
    if abs(l11c) > pr1
        p=0;       
    else 
        L11y1c=cellfun(@(U,V,W,Z) U-l11c*V-l12*W-l13*Z, y1,W1y1,W2y1,W3y1, 'UniformOutput', false);
        f1c=cellfun(@(U,V) U+V, L11y1c,ftemp1, 'UniformOutput', false); 
        f1f1=cellfun(@(U) U'*U, f1c, 'UniformOutput', false);
        f1f2=cellfun(@(U,V) U'*V, f1c,f2, 'UniformOutput', false);
        Fc=F;
        Fc(1,1)=sum(cell2mat(f1f1));
        Fc(1,2)=sum(cell2mat(f1f2));
        Fc(2,1)=Fc(1,2);
        exptermc=0.5*(trace(Fc*psitil));
        L11L22c=cellfun(@(A,B,C,D,E,F,G,H,I,K,L,M,N) ...
            A-(l11c+l21)*B-(l12+l22)*C-(l13+l23)*D ...
            +l11c*l21*E+l11c*l22*F+l11c*l23*G...
            +l12*l21*H+l12*l22*I+l12*l23*K...
            +l13*l21*L+l13*l22*M+l13*l23*N,...
            Ic,W1,W2,W3,W1W1,W1W2,W1W3,W2W1,W2W2,W2W3,W3W1,W3W2,W3W3,...
            'UniformOutput', false); 
        tempD=cellfun(@(U,V) U-V, L11L22c,g1g2Ic, 'UniformOutput', false); 
        detm =  arrayfun(@(i) det(tempD{i,1}),(1:t)','UniformOutput',false);
        detm = log(prod(cell2mat(detm)));  
        cd1c =  detm - exptermc;
        p = exp(cd1c-cd1);
        p = min(1,p);
    end
   
    if p>unifrnd(0,1)
        dl=l11-l11c; 
        l11 = l11c; 
        L11=cellfun(@(U,V) U+dl*V, L11,W1,'UniformOutput',false);
        L11y1=L11y1c;
        accl = accl+1; 
        expterm=exptermc;
        f1=f1c;
        F=Fc;
        cd1=cd1c;
        L11L22=L11L22c;
    end
    
    if accl/m < 0.4
          ccl = ccl/1.1;
    elseif accl/m > 0.6
          ccl = ccl*1.1;
    end  

end


function [l12,L11,L11y1,L11L22,ccl,accl,expterm,f1,F,cd1] = s_p_l12_3W(l12,l11,l13,pr1,L11y1,L11,L11L22,y1,f1,f2,ftemp1,psitil,ccl,accl,expterm,F,m,t,cd1,Ic,g1g2Ic,W1,W2,W3,W1W1,W1W2,W1W3,W2W1,W2W2,W2W3,W3W1,W3W2,W3W3,l21,l22,l23,W1y1,W2y1,W3y1)
     
    % Metropolis Hastings procedure for lambda 12
    
    l12c = l12+ ccl*randn(1,1);
  
    if abs(l12c) > pr1
        p=0;  
    else 
        L11y1c=cellfun(@(U,V,W,Z) U-l11*V-l12c*W-l13*Z, y1,W1y1,W2y1,W3y1, 'UniformOutput', false);
        f1c=cellfun(@(U,V) U+V, L11y1c,ftemp1, 'UniformOutput', false); 
        f1f1=cellfun(@(U) U'*U, f1c, 'UniformOutput', false);
        f1f2=cellfun(@(U,V) U'*V, f1c,f2, 'UniformOutput', false);
        Fc=F;
        Fc(1,1)=sum(cell2mat(f1f1));
        Fc(1,2)=sum(cell2mat(f1f2));
        Fc(2,1)=Fc(1,2);
        exptermc=0.5*(trace(Fc*psitil));
        L11L22c=cellfun(@(A,B,C,D,E,F,G,H,I,K,L,M,N) ...
            A-(l11+l21)*B-(l12c+l22)*C-(l13+l23)*D ...
            +l11*l21*E+l11*l22*F+l11*l23*G...
            +l12c*l21*H+l12c*l22*I+l12c*l23*K...
            +l13*l21*L+l13*l22*M+l13*l23*N,...
            Ic,W1,W2,W3,W1W1,W1W2,W1W3,W2W1,W2W2,W2W3,W3W1,W3W2,W3W3,...
            'UniformOutput', false); 
        tempD=cellfun(@(U,V) U-V, L11L22c,g1g2Ic, 'UniformOutput', false); 
        detm =  arrayfun(@(i) det(tempD{i,1}),(1:t)','UniformOutput',false);
        detm = log(prod(cell2mat(detm)));  
        cd1c =  detm - exptermc;
        p = exp(cd1c-cd1);
        p = min(1,p);
    end
   
    if p>unifrnd(0,1)
        dl=l12-l12c; 
        L11=cellfun(@(U,V) U+dl*V, L11,W2,'UniformOutput',false);
        l12 = l12c;
        L11y1=L11y1c;
        accl = accl+1; 
        expterm=exptermc;
        f1=f1c;
        F=Fc;
        cd1=cd1c;
        L11L22=L11L22c;
    end
    
    if accl/m < 0.4
          ccl = ccl/1.1;
    elseif accl/m > 0.6
          ccl = ccl*1.1;
    end  
 
end


function [l13,L11,L11y1,L11L22,ccl,accl,expterm,f1,F,cd1] = s_p_l13_3W(l13,l11,l12,pr1,L11y1,L11,L11L22,y1,f1,f2,ftemp1,psitil,ccl,accl,expterm,F,m,t,cd1,Ic,g1g2Ic,W1,W2,W3,W1W1,W1W2,W1W3,W2W1,W2W2,W2W3,W3W1,W3W2,W3W3,l21,l22,l23,W1y1,W2y1,W3y1)
     
    % Metropolis Hastings procedure for lambda 13
    
    l13c = l13+ ccl*randn(1,1);
  
    if abs(l13c) > pr1
        p=0;  
    else 
        L11y1c=cellfun(@(U,V,W,Z) U-l11*V-l12*W-l13c*Z, y1,W1y1,W2y1,W3y1, 'UniformOutput', false);
        f1c=cellfun(@(U,V) U+V, L11y1c,ftemp1, 'UniformOutput', false); 
        f1f1=cellfun(@(U) U'*U, f1c, 'UniformOutput', false);
        f1f2=cellfun(@(U,V) U'*V, f1c,f2, 'UniformOutput', false);
        Fc=F;
        Fc(1,1)=sum(cell2mat(f1f1));
        Fc(1,2)=sum(cell2mat(f1f2));
        Fc(2,1)=Fc(1,2);
        exptermc=0.5*(trace(Fc*psitil));
        L11L22c=cellfun(@(A,B,C,D,E,F,G,H,I,K,L,M,N) ...
            A-(l11+l21)*B-(l12+l22)*C-(l13c+l23)*D ...
            +l11*l21*E+l11*l22*F+l11*l23*G...
            +l12*l21*H+l12*l22*I+l12*l23*K...
            +l13c*l21*L+l13c*l22*M+l13c*l23*N,...
            Ic,W1,W2,W3,W1W1,W1W2,W1W3,W2W1,W2W2,W2W3,W3W1,W3W2,W3W3,...
            'UniformOutput', false);  
        tempD=cellfun(@(U,V) U-V, L11L22c,g1g2Ic, 'UniformOutput', false); 
        detm =  arrayfun(@(i) det(tempD{i,1}),(1:t)','UniformOutput',false);
        detm = log(prod(cell2mat(detm)));  
        cd1c =  detm - exptermc;
        p = exp(cd1c-cd1);
        p = min(1,p);
    end
   
    if p>unifrnd(0,1)
        dl=l13-l13c; 
        L11=cellfun(@(U,V) U+dl*V, L11,W3,'UniformOutput',false);
        l13 = l13c;
        L11y1=L11y1c;
        accl = accl+1; 
        expterm=exptermc;
        f1=f1c;
        F=Fc;
        cd1=cd1c;
        L11L22=L11L22c;
    end
    
    if accl/m < 0.4
          ccl = ccl/1.1;
    elseif accl/m > 0.6
          ccl = ccl*1.1;
    end  
 
end


function [l21,L22,L22y2,L11L22,ccl,accl,expterm,f2,F,cd1] = s_p_l21_3W(l21,l22,l23,pr1,L22y2,L22,L11L22,y2,f1,f2,ftemp2,psitil,ccl,accl,expterm,F,m,t,cd1,Ic,g1g2Ic,W1,W2,W3,W1W1,W1W2,W1W3,W2W1,W2W2,W2W3,W3W1,W3W2,W3W3,l11,l12,l13,W1y2,W2y2,W3y2)

    % Metropolis Hastings procedure for lambda 21
    
    l21c = l21+ ccl*randn(1,1);
    
    if abs(l21c) > pr1
        p=0;    
    else
        L22y2c=cellfun(@(U,V,W,Z) U-l21c*V-l22*W-l23*Z, y2,W1y2,W2y2,W3y2, 'UniformOutput', false);
        f2c=cellfun(@(U,V) U+V, L22y2c,ftemp2, 'UniformOutput', false); 
        f1f2=cellfun(@(U,V) U'*V, f1,f2c, 'UniformOutput', false);
        f2f2=cellfun(@(U) U'*U, f2c, 'UniformOutput', false); 
        Fc=F;
        Fc(1,2)=sum(cell2mat(f1f2));
        Fc(2,1)=Fc(1,2);
        Fc(2,2)=sum(cell2mat(f2f2));
        exptermc=0.5*(trace(Fc*psitil));
        L11L22c=cellfun(@(A,B,C,D,E,F,G,H,I,K,L,M,N) ...
            A-(l11+l21c)*B-(l12+l22)*C-(l13+l23)*D ...
            +l11*l21c*E+l11*l22*F+l11*l23*G...
            +l12*l21c*H+l12*l22*I+l12*l23*K...
            +l13*l21c*L+l13*l22*M+l13*l23*N,...
            Ic,W1,W2,W3,W1W1,W1W2,W1W3,W2W1,W2W2,W2W3,W3W1,W3W2,W3W3,...
            'UniformOutput', false); 
        tempD=cellfun(@(U,V) U-V, L11L22c,g1g2Ic, 'UniformOutput', false); 
        detm =  arrayfun(@(i) det(tempD{i,1}),(1:t)','UniformOutput',false);
        detm =  log(prod(cell2mat(detm)));  
        cd1c =  detm - exptermc;
        p = exp(cd1c-cd1);
        p = min(1,p);
    end
   
    if p>unifrnd(0,1)
        dl=l21-l21c; 
        l21 = l21c; 
        L22=cellfun(@(U,V) U+dl*V, L22,W1,'UniformOutput',false);
        L22y2=L22y2c;
        L11L22=L11L22c;
        accl = accl+1; 
        expterm=exptermc;
        f2=f2c;
        F=Fc;
        cd1=cd1c;
    end

    if accl/m < 0.4
          ccl = ccl/1.1;
    elseif accl/m > 0.6
          ccl = ccl*1.1;
    end  
    
end


function [l22,L22,L22y2,L11L22,ccl,accl,expterm,f2,F,cd1] = s_p_l22_3W(l22,l21,l23,pr1,L22y2,L22,L11L22,y2,f1,f2,ftemp2,psitil,ccl,accl,expterm,F,m,t,cd1,Ic,g1g2Ic,W1,W2,W3,W1W1,W1W2,W1W3,W2W1,W2W2,W2W3,W3W1,W3W2,W3W3,l11,l12,l13,W1y2,W2y2,W3y2)
    
    % Metropolis Hastings procedure for lambda 22
    
    l22c = l22+ ccl*randn(1,1);
    
    if abs(l22c) > pr1
        p=0;    
    else
        L22y2c=cellfun(@(U,V,W,Z) U-l21*V-l22c*W-l23*Z, y2,W1y2,W2y2,W3y2, 'UniformOutput', false);
        f2c=cellfun(@(U,V) U+V, L22y2c,ftemp2, 'UniformOutput', false); 
        f1f2=cellfun(@(U,V) U'*V, f1,f2c, 'UniformOutput', false);
        f2f2=cellfun(@(U) U'*U, f2c, 'UniformOutput', false); 
        Fc=F;
        Fc(1,2)=sum(cell2mat(f1f2));
        Fc(2,1)=Fc(1,2);
        Fc(2,2)=sum(cell2mat(f2f2));
        exptermc=0.5*(trace(Fc*psitil));
        L11L22c=cellfun(@(A,B,C,D,E,F,G,H,I,K,L,M,N) ...
            A-(l11+l21)*B-(l12+l22c)*C-(l13+l23)*D ...
            +l11*l21*E+l11*l22c*F+l11*l23*G...
            +l12*l21*H+l12*l22c*I+l12*l23*K...
            +l13*l21*L+l13*l22c*M+l13*l23*N,...
            Ic,W1,W2,W3,W1W1,W1W2,W1W3,W2W1,W2W2,W2W3,W3W1,W3W2,W3W3,...
            'UniformOutput', false); 
        tempD=cellfun(@(U,V) U-V, L11L22c,g1g2Ic, 'UniformOutput', false); 
        detm =  arrayfun(@(i) det(tempD{i,1}),(1:t)','UniformOutput',false);
        detm =  log(prod(cell2mat(detm)));  
        cd1c =  detm - exptermc;
        p = exp(cd1c-cd1);
        p = min(1,p);
    end
   
    if p>unifrnd(0,1)
        dl=l22-l22c; 
        l22 = l22c; 
        L22=cellfun(@(U,V) U+dl*V, L22,W2,'UniformOutput',false);
        L22y2=L22y2c;
        L11L22=L11L22c;
        accl = accl+1; 
        expterm=exptermc;
        f2=f2c;
        F=Fc;
        cd1=cd1c;
    end

    if accl/m < 0.4
          ccl = ccl/1.1;
    elseif accl/m > 0.6
          ccl = ccl*1.1;
    end  
    
end


function [l23,L22,L22y2,L11L22,ccl,accl,expterm,f2,F,cd1] = s_p_l23_3W(l23,l21,l22,pr1,L22y2,L22,L11L22,y2,f1,f2,ftemp2,psitil,ccl,accl,expterm,F,m,t,cd1,Ic,g1g2Ic,W1,W2,W3,W1W1,W1W2,W1W3,W2W1,W2W2,W2W3,W3W1,W3W2,W3W3,l11,l12,l13,W1y2,W2y2,W3y2)

    % Metropolis Hastings procedure for lambda 23
    
    l23c = l23+ ccl*randn(1,1);
    
    if abs(l23c) > pr1  
        p=0;       
    else
        L22y2c=cellfun(@(U,V,W,Z) U-l21*V-l22*W-l23c*Z, y2,W1y2,W2y2,W3y2, 'UniformOutput', false);
        f2c=cellfun(@(U,V) U+V, L22y2c,ftemp2, 'UniformOutput', false); 
        f1f2=cellfun(@(U,V) U'*V, f1,f2c, 'UniformOutput', false);
        f2f2=cellfun(@(U) U'*U, f2c, 'UniformOutput', false); 
        Fc=F;
        Fc(1,2)=sum(cell2mat(f1f2));
        Fc(2,1)=Fc(1,2);
        Fc(2,2)=sum(cell2mat(f2f2));
        exptermc=0.5*(trace(Fc*psitil));
        L11L22c=cellfun(@(A,B,C,D,E,F,G,H,I,K,L,M,N) ...
            A-(l11+l21)*B-(l12+l22)*C-(l13+l23c)*D ...
            +l11*l21*E+l11*l22*F+l11*l23c*G...
            +l12*l21*H+l12*l22*I+l12*l23c*K...
            +l13*l21*L+l13*l22*M+l13*l23c*N,...
            Ic,W1,W2,W3,W1W1,W1W2,W1W3,W2W1,W2W2,W2W3,W3W1,W3W2,W3W3,...
            'UniformOutput', false); 
        tempD=cellfun(@(U,V) U-V, L11L22c,g1g2Ic, 'UniformOutput', false); 
        detm =  arrayfun(@(i) det(tempD{i,1}),(1:t)','UniformOutput',false);
        detm =  log(prod(cell2mat(detm)));  
        cd1c =  detm - exptermc;
        p = exp(cd1c-cd1);
        p = min(1,p);
    end
   
    if p>unifrnd(0,1)
        dl=l23-l23c; 
        l23 = l23c; 
        L22=cellfun(@(U,V) U+dl*V, L22,W3,'UniformOutput',false);
        L22y2=L22y2c;
        L11L22=L11L22c;
        accl = accl+1; 
        expterm=exptermc;
        f2=f2c;
        F=Fc;
        cd1=cd1c;
    end

    if accl/m < 0.4
          ccl = ccl/1.1;
    elseif accl/m > 0.6
          ccl = ccl*1.1;
    end  

end
