function [lmeanvarPvalue lindPvalue lbothPvalue] = berkowitz(y)
% modified version of Berkowitz (2001) codes, since then fmins has been
% substituted by fminsearch in matlab

    capt = length(y);
    for ll = 1:capt
        y(ll)=min(y(ll),.999); % get rid of corner values
        y(ll)=max(y(ll),.001);
        z(ll)=norminv(y(ll));
    end

    z=z(:);

    %%% BASIC LR METHOD
    oparam = [.02 .02];
    options.MaxIter = 1000;
    [minp, minval] = fminsearch('carloss',oparam,options,z,capt);
    fmin = -minval;
    fmeanvar = -lossvar([0,minp(2),1] , z,capt);
    findep  = -carloss([minp(1),0], z, capt);
    fboth   = -lossvar([0,0,1], z, capt);

    lmeanvar    =2*(fmin-fmeanvar);
    lind        =2*(fmin-findep);
    lboth       =2*(fmin-fboth);

    lmeanvarPvalue =  1-cdf('chi2',lmeanvar,2);
    lindPvalue = 1-cdf('chi2',lind,1);
    lbothPvalue = 1-cdf('chi2',lboth,3);
end


function [f,g] = lossvar(param,y,capt)

    mu=param(1);
    rho=param(2);
    sig=param(3);
    ly=lx(y);
    etahat=y-mu-rho.*ly;
    etahat=etahat(2:capt,:);


    f=-.5*log(2*pi)-.5*log(sig/(1-rho^2)) ...
      -(y(1)-(mu/(1-rho))^2)/(2*sig/(1-rho^2)) ...
      -(capt-1)/2*log(2*pi) ...
      -(capt-1)/2*log(sig)-sum((etahat.^2)./(2*sig));

    f=-f;
end