function ordered_probit_empirical(dataset,ntest,indx_train,indx_test)       
%%%%%%%%%%%%%%%%%%%%%%  Part I: estimate the coefficients and the probability forecasts of all sub-models  %%%%%%%%%%%%%%% 
warning off;
samptest  = dataset(indx_test,:);%test sample
samptrain = dataset(indx_train,:);%training sample
[nnew,vartest] = size(samptest);
Param0    = [-0.02;-0.01;0.01;0.02;[0.001;0.0100;0.05;-0.002;-0.001;-0.001;0.001;0.01;0.001]];%The initial value of coefficients that is input to the maximum likelihood function.
options   = optimset('display','off','MaxIter',1000,'MaxFunEvals',3500,'TolX',10^-5,'TolFun',10^-5);
zerovalue = 0.0001;
sall      = length(Param0);%sall is the number of all variables
snull     = 4;%snull is the number of mandatory variables
pfocus    = snull;
puncert   = (sall-snull);%puncert is the number of alternative variables
I         = eye(puncert);
[s,S]     = Choose_variable(1,sall-snull+1);%Choose_variable is a function used to choose variables for each sub-model.
temp      = S-(S>0);
S0        = temp(:,2:sall-snull+1);
M         = length(s);%M is the number of all alternative models.
w0        = ones(M,1)/M;
tic;
worry=0;
ycode     = samptrain(:,1:5);
ytrue     = samptrain(:,6);
X         = samptrain(:,7:15);
Xnew      = samptest(:,7:15);
ynew      = samptest(:,6);
Xtrain    = X;
n         = length(ytrue);
logn      = log(n);
nroot     = n^0.5;
X         = [ones(n,1),X];
betahat   = zeros(sall,M);
    for m=1:M
       p            = s(m,1);
       a            = S(m,1:p);
       Xm           = X(:,a);
       param0       = Param0(1:snull-1,:);
       param0       = [param0;Param0(snull-1+a,:)];
       [paramhat hessian]=mle_ordered_probit(param0,ycode,Xm);%estimate coefficients by maximum likelihood estimation method
       [lc(m,1)]    = con_log_l_ordered_probit(paramhat,ycode,Xm);% conditional log-likelihood function of ordinal data
       LH(m,1)      = lc(m,1);
       AIC(m,1)     = 2*LH(m,1)+2*(3+p);%AIC value of the m-th alternative model
       BIC(m,1)     = 2*LH(m,1)+logn*(3+p);%BIC value of the m-th alternative model
       betahat([1;2;3;3+a'],m) = paramhat;%coefficient estimations of the m-th alternative model
       alpha1hat    = betahat(1,m);
       alpha2hat    = betahat(2,m);
       alpha3hat    = betahat(3,m);
       alpha4hat    = betahat(4,m);
       Xtrainbetahat= Xtrain*betahat(5:sall,m);
       X1train      = Xtrainbetahat+alpha1hat;
       X2train      = Xtrainbetahat+alpha2hat;
       X3train      = Xtrainbetahat+alpha3hat;
       X4train      = Xtrainbetahat+alpha4hat;
       F1train      = normcdf(X1train);
       F2train      = normcdf(X2train);
       F3train      = normcdf(X3train);
       F4train      = normcdf(X4train);       
       Xbetahat     = Xnew*betahat(5:sall,m);
       X1           = Xbetahat+alpha1hat;
       X2           = Xbetahat+alpha2hat;
       X3           = Xbetahat+alpha3hat;
       X4           = Xbetahat+alpha4hat;      
       F1           = normcdf(X1);
       F2           = normcdf(X2);
       F3           = normcdf(X3);
       F4           = normcdf(X4);
       P(:,m,5)     = (1-F4);% the probability forcasts of Choice 5 from the m-th alternative model for all observations in the test sample
       P(:,m,4)     = (F4-F3);% the probability forcasts of Choice 4 from the m-th alternative model for all observations in the test sample
       P(:,m,3)     = (F3-F2);% the probability forcasts of Choice 3 from the m-th alternative model for all observations in the test sample
       P(:,m,2)     = (F2-F1);% the probability forcasts of Choice 2 from the m-th alternative model for all observations in the test sample
       P(:,m,1)     = F1;% the probability forcasts of Choice 1 from the m-th alternative model for all observations in the test sample
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Part I ends here %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
    %%
    [screening,Nscreening]    = topfivescreening(BIC,S);%topfivescreening function is used to choose the best 5 five sub-models according to BIC
    ns                = length(Nscreening);
    %% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Part II£ºmodel averaging probability forecasts by using LZWZ, A-opt and S-FIC %%%%%%%%%%%%%%%%%%%%%    
    J                = hessian/n;
    J00              = J(1:pfocus,1:pfocus);
    J01              = J(1:pfocus,pfocus+1:sall);
    J10              = J01';
    J11              = J(pfocus+1:sall,pfocus+1:sall);
    invJ00           = eye(pfocus)/J00;
    J10invJ00        = J10*invJ00;
    Kinv             = J11-J10invJ00*J01;
    delta            = nroot*paramhat(pfocus+1:sall,1);
    Kinvdelta        = Kinv*delta;
    for m=2:M
        p           = s(m,1)-1;
        a           = S0(m,1:p);
        pii         = I(a,:);
        Km(:,:,m)   = pii'/(pii*Kinv*pii')*pii;
        tt(:,m)     = delta-Km(:,:,m)*Kinvdelta;
    end
    tt(:,1)         = delta;
    %%
    %omega is getted by letting the partial derivatives evaluated at the null point
    X1              = alpha1hat+zeros(nnew,1);
    X2              = alpha2hat+zeros(nnew,1);
    X3              = alpha3hat+zeros(nnew,1);
    X4              = alpha4hat+zeros(nnew,1);
    normpdfnullX1       = normpdf(X1);
    normpdfnullX2       = normpdf(X2);
    normpdfnullX3       = normpdf(X3);
    normpdfnullX4       = normpdf(X4);  
    j            = 1;
    while j<=nnew
        for fi=1:5
            if fi==1
                diffmutheta = [normpdfnullX1(j,1);0;0;0];
                diffmugamma = normpdfnullX1(j,1)*Xnew(j,:)';
            end
            if fi==2
                diffmutheta = [-normpdfnullX1(j,1);normpdfnullX2(j,1);0;0];
                diffmugamma = (normpdfnullX2(j,1)-normpdfnullX1(j,1))*Xnew(j,:)';
            end
            if fi==3
                diffmutheta = [0;-normpdfnullX2(j,1);normpdfnullX3(j,1);0];
                diffmugamma = (normpdfnullX3(j,1)-normpdfnullX2(j,1))*Xnew(j,:)';                    
            end
            if fi==4
                diffmutheta = [0;0;-normpdfnullX3(j,1);normpdfnullX4(j,1)];
                diffmugamma = (normpdfnullX4(j,1)-normpdfnullX3(j,1))*Xnew(j,:)';                    
            end
            if fi==5
                diffmutheta = [0;0;0;-normpdfnullX4(j,1)];
                diffmugamma = -normpdfnullX4(j,1)*Xnew(j,:)';
            end
            omega       = J10invJ00*diffmutheta-diffmugamma;
                  
            T           = zeros(M,M);
            for m=1:M
                Kmomega(:,m)     = Km(:,:,m) * omega;
                for m1=1:m
                    T(m,m1)      = Kmomega(:,m)'*Kinv* Kmomega(:,m1);
                end
                t(m,1)       = omega'*tt(:,m);
            end
            T                = T+T'-diag(diag(T));
            TT               = T+t*t';
            w                = quadprog(TT,zeros(1,M),zeros(1,M),0,ones(1,M),1,zeros(M,1),ones(M,1),w0,options);%w is the weight vector derived from A-opt method
            Pma(j,fi)        = P(j,:,fi)*w;%Pma(j,fi) is the model averaging probability forecasts of Choice fi for the j-th observation in the test sample by using A-opt method
            %%
            omegaKmomega    = Kmomega'*omega;
            TT_v2           = t*t';
            w_v2            = quadprog(TT_v2,omegaKmomega,zeros(1,M),0,ones(1,M),1,zeros(M,1),ones(M,1),w0,options);%w_v2 is the weight vector derived from LZWZ method     
            Pma_v2(j,fi)    = P(j,:,fi)*w_v2;%Pma_v2(j,fi)is the model averaging probability forecasts of Choice fi for the j-th observation in the test sample by using LZWZ method
            %%
            FIC             = (diag(TT_v2)+2*omegaKmomega)/omegaKmomega(M,1);%FIC is the vector includes the FIC value of each alternative models
            w               = zeros(M,1);
            wFIC            = w;
            [temp,xic]      = min(FIC);
            [~,ms_FIC]      = min(FIC);
            wFIC(xic,1)     = 1;
            FIC             = FIC-min(FIC);
            sFIC            = exp(-FIC/2)/sum(exp(-FIC/2));%sFIC is the weight vector derived from S-FIC without screening method
            sFICscreening   = exp(-FIC(Nscreening',1)/2)/sum(exp(-FIC(Nscreening',1)/2));%sFICscreening is the weight vector derived from S-FIC with screening method
            PFIC(j,fi)      = P(j,:,fi)*wFIC;%the model selection probability forecasts of Choice fi for the j-th observation in the test sample by using FIC as model selection criterion
            PsFIC(j,fi)           = P(j,:,fi)*sFIC;%the model averaging probability forecasts of Choice fi for the j-th observation in the test sample by using S-FIC without screening method
            PsFICscreening(j,fi)  = P(j,Nscreening,fi)*sFICscreening;%the model averaging probability forecasts of Choice fi for the j-th observation in the test sample by using S-FIC with screening method
            %% screening
            TTscreening       = TT(Nscreening',Nscreening);
            TT_v2screening    = TT_v2(Nscreening',Nscreening);
            wscreening        = quadprog(TTscreening,zeros(1,ns),zeros(1,ns),0,ones(1,ns),1,zeros(ns,1),ones(ns,1),ones(ns,1)/ns,options);%wscreening is the weight vector derived from A-opt with screening method
            w_v2screening     = quadprog(TT_v2screening,omegaKmomega(Nscreening,1),zeros(1,ns),0,ones(1,ns),1,zeros(ns,1),ones(ns,1),ones(ns,1)/ns,options);%wscreening is the weight vector derived from LZWZ with screening method
            Pmascreening(j,fi)        = P(j,Nscreening,fi)*wscreening;%the model averaging probability forecasts of Choice fi for the j-th observation in the test sample by using A-opt with screening method
            Pma_v2screening(j,fi)     = P(j,Nscreening,fi)*w_v2screening;%the model averaging probability forecasts of Choice fi for the j-th observation in the test sample by using LZWZ with screening method
            if max(abs(w-w0))<zerovalue  %|| max(abs(wtrue-w0))<zerovalue|| max(abs(wscreening-ones(ns,1)/(ns)))<zerovalue || max(abs(wtruescreening-ones(ns,1)/(ns)))<zerovalue
                worry=1;
            end
        end
        if worry==0
            j=j+1;
        else
            j=nnew+1;
        end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Part II ends here %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Part III: model averaging probability forecasts by using JMA method %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    A2           = zeros(M,M);
    Xc           = Xtrain;
    cbetahat     = zeros(sall,M);
    for i=1:n
       cX=[X(1:i-1,:);X(i+1:n,:)];
       cycode=[ycode(1:i-1,:);ycode(i+1:n,:)];
       for m=1:M      
           p            = s(m,1);
           a            = S(m,1:p);
           cXm          = cX(:,a);
           param0       = Param0(1:snull-1,:);
           param0       = [param0;Param0(snull-1+a,:)];
           [cparamhat hessian]  = mle_ordered_probit(param0,cycode,cXm);%cparamhat is the coefficient estimators for the m-th alternative model by deleting the i-th observation from the training sample
           cbetahat([1;2;3;3+a'],m) = cparamhat;
           alpha1hat    = cbetahat(1,m);
           alpha2hat    = cbetahat(2,m);
           alpha3hat    = cbetahat(3,m);
           alpha4hat    = cbetahat(4,m);
           Xbetahat     = Xc(i,:)*cbetahat(5:sall,m);
           X1           = Xbetahat+alpha1hat;
           X2           = Xbetahat+alpha2hat;
           X3           = Xbetahat+alpha3hat;
           X4           = Xbetahat+alpha4hat;        
           F1           = normcdf(X1);
           F2           = normcdf(X2);
           F3           = normcdf(X3);
           F4           = normcdf(X4);
           PP(m,5)      = (1-F4);%the probability forecast of Choice 5 for the i-th observation based on the m-th sub-model which is estimated using the training samplewith the i-th observation deleted.
           PP(m,4)      = (F4-F3);%the probability forecast of Choice 4 for the i-th observation based on the m-th sub-model which is estimated using the training samplewith the i-th observation deleted.
           PP(m,3)      = (F3-F2);%the probability forecast of Choice 3 for the i-th observation based on the m-th sub-model which is estimated using the training samplewith the i-th observation deleted.
           PP(m,2)      = (F2-F1);%the probability forecast of Choice 2 for the i-th observation based on the m-th sub-model which is estimated using the training samplewith the i-th observation deleted.
           PP(m,1)      = F1;%the probability forecast of Choice 1 for the i-th observation based on the m-th sub-model which is estimated using the training samplewith the i-th observation deleted.
           Pfinal(m,1)  = PP(m,:)*ycode(i,:)';
       end  
       B2=(PP(:,1)-ycode(i,1))*(PP(:,1)-ycode(i,1))'+(PP(:,2)-ycode(i,2))*(PP(:,2)-ycode(i,2))'...
                +(PP(:,3)-ycode(i,3))*(PP(:,3)-ycode(i,3))'+(PP(:,4)-ycode(i,4))*(PP(:,4)-ycode(i,4))'+(PP(:,5)-ycode(i,5))*(PP(:,5)-ycode(i,5))';
       A2=A2+B2;
    end
    wcross2=quadprog(A2,zeros(M,1),zeros(1,M),0,ones(1,M),1,zeros(M,1),ones(M,1),w0,options);%wcross2 is the weight vector derived from the JMA without screening method.
    A2screening=A2(Nscreening',Nscreening);
    wcrossscreening2=quadprog(A2screening,zeros(ns,1),zeros(1,ns),0,ones(1,ns),1,zeros(ns,1),ones(ns,1),ones(ns,1)/ns,options);%wcrossscreening2 is the weight vector derived from the JMA with screening method.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Part III ends here %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Part IV: model averaging probability forecasts by using S-AIC£¬ S-BIC, EW method%%%%%%%%%%%%%%%%%%%%%% 
w           = zeros(M,1);
wAIC        = w;
[temp,xic]  = min(AIC);
[~,ms_AIC]  = min(AIC);
wAIC(xic,1) = 1;
AIC         = AIC-min(AIC);
sAIC        = exp(-AIC/2)/sum(exp(-AIC/2));%sAIC is the weight vector derived from the S-AIC without screening method.
sAICscreening             = exp(-AIC(Nscreening',1)/2)/sum(exp(-AIC(Nscreening',1)/2));%sAICscreening is the weight vector derived from the S-AIC with screening method.

wBIC        = w;
[temp,xic]  = min(BIC);
[~,ms_BIC]  = min(BIC);
wBIC(xic,1) = 1;
BIC         = BIC-min(BIC);
sBIC        = exp(-BIC/2)/sum(exp(-BIC/2));%sBIC is the weight vector derived from the S-BIC without screening method.
sBICscreening             = exp(-BIC(Nscreening',1)/2)/sum(exp(-BIC(Nscreening',1)/2));%sBICscreening is the weight vector derived from the S-BIC with screening method.
    for fi=1:5
        PAIC(:,fi)                   = P(:,:,fi)*wAIC;%the probability forecast of Choice fi by using AIC method
        PsAIC(:,fi)                  = P(:,:,fi)*sAIC;%the probability forecast of Choice fi by using sAIC without screening method
        PsAICscreening(:,fi)         = P(:,Nscreening,fi)*sAICscreening;%the probability forecast of Choice fi by using sAIC with screening method
        PBIC(:,fi)                   = P(:,:,fi)*wBIC;%the probability forecast of Choice fi by using BIC method
        PsBIC(:,fi)                  = P(:,:,fi)*sBIC;%the probability forecast of Choice fi by using sBIC without screening method
        PsBICscreening(:,fi)         = P(:,Nscreening,fi)*sBICscreening;%the probability forecast of Choice fi by using sBIC with screening method
        Pequal(:,fi)                 = P(:,:,fi)*ones(M,1)/M;%the probability forecast of Choice fi by using EW without screening method
        Pequalscreening(:,fi)        = P(:,Nscreening,fi)*ones(ns,1)/ns;%the probability forecast of Choice fi by using EW with screening method
        Pcross2(:,fi)                 = P(:,:,fi)*wcross2;%the probability forecast of Choice fi by using JMA without screening method
        Pcrossscreening2(:,fi)        = P(:,Nscreening,fi)*wcrossscreening2;%the probability forecast of Choice fi by using JMA with screening method
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Part IV ends here %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Part V: Output from all model selection and model averaging method%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    PALL(:,:,1)  =  PAIC;                 PALL(:,:,2)  =  PBIC;               PALL(:,:,3)  =  PFIC;     
    PALL(:,:,4)  =  PsAIC;                PALL(:,:,5)  =  PsBIC;              PALL(:,:,6)  =  PsFIC;            
    PALL(:,:,7)  =  Pequal;               PALL(:,:,8) =  Pma_v2;              PALL(:,:,9)  =  Pma;                          
    PALL(:,:,10) =  PsAICscreening;       PALL(:,:,11) =  PsBICscreening;     PALL(:,:,12) =  PsFICscreening;
    PALL(:,:,13) =  Pequalscreening;      PALL(:,:,14) =  Pma_v2screening;    PALL(:,:,15) =  Pmascreening;
    PALL(:,:,16) =  Pcross2;              PALL(:,:,17) =  Pcrossscreening2;
    %
    for pa=1:17
        [Result(1,pa)] = fun2(ynew,PALL(:,:,pa));%fun2 is used to calculate the hit rate of all model selection and model averaging method
    end
Resultall = Result;
xlswrite(['empirical_ordered_probit_ntest=',num2str(ntest)],Resultall);
ModelSlected = [ms_AIC ms_BIC ms_FIC];%the models that AIC BIC and FIC select seperately
xlswrite(['ModelSlected_ntest=',num2str(ntest)],ModelSlected);
betahat_MAsBIC             =  betahat*sBIC;%the model generated by combining the coefficient estimators of sub-models using S-BIC without screening method.
betahat_MAsAIC             =  betahat*sAIC;%the model generated by combining the coefficient estimators of sub-models using S-AIC without screening method.
betahat_MAEW               =  betahat*ones(M,1)/M;%the model generated by combining the coefficient estimators of sub-models using EW without screening method.
betahat_MAJMA              =  betahat*wcross2;%the model generated by combining the coefficient estimators of sub-models using JMA without screening method.
betahat_MAsBIC_screening       =  betahat(:,Nscreening)*sBICscreening;%the model generated by combining the coefficient estimators of sub-models using S-BIC with screening method.
betahat_MAsAIC_screening       =  betahat(:,Nscreening)*sAICscreening;%the model generated by combining the coefficient estimators of sub-models using S-AIC with screening method.
betahat_MAEW_screening         =  betahat(:,Nscreening)*ones(5,1)/5;%the model generated by combining the coefficient estimators of sub-models using EW with screening method.
betahat_MAJMA_screening        =  betahat(:,Nscreening)*wcrossscreening2;%the model generated by combining the coefficient estimators of sub-models using JMA with screening method.
betahat_all               =  [betahat,betahat_MAsBIC,betahat_MAsAIC,betahat_MAEW,betahat_MAJMA,betahat_MAsBIC_screening,betahat_MAsAIC_screening,betahat_MAEW_screening,betahat_MAJMA_screening];
xlswrite(['coefficient_estimation=',num2str(ntest),'.xlsx'],betahat_all);
toc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Part V ends here %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [correct]  = fun2(ynew,Pw)
    [temp,yhat]    = max(Pw');
    yhat           = yhat';
    correct        = mean(yhat==ynew);
end

function [paramhat hessian]=mle_ordered_probit(param0,ycode,X)% The maximum likelihood estimation function to estimate the coefficients of sub-models.
d1=1;
options = optimset('Display','off'); 
err=0.0001;
while d1>err
    [paramhat, fvalue, EXITFLAG,OUTPUT,hessian]=fsolve(@(param) Ob_Funwithoutgrad(param,X,ycode),param0,options);
    max(abs(fvalue));
    d1=norm(paramhat-param0);
    param0=paramhat;
end
function [Ob_Funwithoutgrad]=Ob_Funwithoutgrad(param,X,ycode)
pp       = length(param);
beta     = param(5:pp,1);
alpha1   = param(1,1);
alpha2   = param(2,1);
alpha3   = param(3,1);
alpha4   = param(4,1);
Xbeta    = X(:,2:pp-3)*beta;
normX1   = Xbeta+alpha1;
normX2   = Xbeta+alpha2;
normX3   = Xbeta+alpha3;
normX4   = Xbeta+alpha4;
normpdfX1= normpdf(normX1);
normpdfX2= normpdf(normX2);
normpdfX3= normpdf(normX3);
normpdfX4= normpdf(normX4);
F1       = normcdf(normX1);
F2       = normcdf(normX2);
F3       = normcdf(normX3);
F4       = normcdf(normX4);
Ptrain   = [F1,F2-F1,F3-F2,F4-F3,1-F4];
diffalpha1= 0;
diffalpha2= 0;
diffalpha3= 0;
diffalpha4= 0;
diffbeta = zeros(pp-4,1);
n        = length(normX1);
for l=1:n
    diffalpha1=diffalpha1+ycode(l,:)./Ptrain(l,:)*[normpdfX1(l,:);-normpdfX1(l,:);0;0;0];
    diffalpha2=diffalpha2+ycode(l,:)./Ptrain(l,:)*[0;normpdfX2(l,:);-normpdfX2(l,:);0;0];
    diffalpha3=diffalpha3+ycode(l,:)./Ptrain(l,:)*[0;0;normpdfX3(l,:);-normpdfX3(l,:);0];
    diffalpha4=diffalpha4+ycode(l,:)./Ptrain(l,:)*[0;0;0;normpdfX4(l,:);-normpdfX4(l,:)];
    diffbeta=diffbeta+ycode(l,:)./Ptrain(l,:)*[normpdfX1(l,:);normpdfX2(l,:)-normpdfX1(l,:);normpdfX3(l,:)-normpdfX2(l,:);normpdfX4(l,:)-normpdfX3(l,:);-normpdfX4(l,:)]*X(l,2:pp-3)';
end
Ob_Funwithoutgrad=[diffalpha1;diffalpha2;diffalpha3;diffalpha4;diffbeta];
end
end

function [lc]=con_log_l_ordered_probit(paramhat,ycode,X)% conditional log-likelihood function of ordinal data
pp     = length(paramhat);
beta   = paramhat(5:pp,1);
alpha1 = paramhat(1,1);
alpha2 = paramhat(2,1);
alpha3 = paramhat(3,1);
alpha4 = paramhat(4,1);
Xbeta  = X(:,2:pp-3)*beta;
normX1 = Xbeta+alpha1;
normX2 = Xbeta+alpha2;
normX3 = Xbeta+alpha3;
normX4 = Xbeta+alpha4;
F1     = normcdf(normX1);
F2     = normcdf(normX2);
F3     = normcdf(normX3);
F4     = normcdf(normX4);
logP5  = log(1-F4);
logP4  = log(F4-F3);
logP3  = log(F3-F2);
logP2  = log(F2-F1);
logP1  = log(F1);
lc     = -sum(sum([logP1 logP2 logP3 logP4 logP5].*ycode));
end
function [s,S]=Choose_variable(k,m)
q=m-k;
M=2^q;
pp=zeros(q+1,1);
for i=1:q+1
    pp(i,1)=nchoosek(q,i-1);
end
s=zeros(M,1);
s(1,1)=k;
for i=2:q+1
    temp=sum(s>0);
    s(temp+1:temp+pp(i,1),1)=(i+k-1)*ones(pp(i,1),1);
end
S=zeros(M,m);
for i=1:M
    S(i,1:k)=[1:k];
end
tS=1;
for i=1:q
    Temp=nchoosek([k+1:m],i);
    S(tS+1:tS+size(Temp,1),k+1:k+i)=Temp;
    tS=tS+size(Temp,1);
end
end
function [ICscreening,NICscreening] = topfivescreening(IC,atemp)

    [ICscreening(1,1),NICscreening(1,1)]=min(IC);
    for i=2:5
        ICi = IC.*(IC>(ICscreening(i-1,1)+eps));
        ICi = ICi+(ICi==0)*eps^(-1);
        [ICscreening(i,1),NICscreening(i,1)]=min(ICi);
    end
end
end