%/*** Procedure to run tests of predictability ***/

function [return1, return2] = predtest(ret,pred,maxp,DFbelt,DFGLSbelt, tconf, tquantL, tquantR, Qconf)

% local zcrit  p  n T r, x, Dx, Lx, Z, i, 
%     coef, v, Sv2, varmat, SErho, e, psi, df, 
%     rhogls, xbar, inter, mu, Zbar, resid, dfgls, 
%     conf, u, beta, Sigma, Kurt, delta, norm, omega2, 
%     normbeta, SEbeta, tstat, tbeta, cL, cU, 
%     critL, critR, Qbeta, rhoL, rhoU, correct, rbar;

% global DFbelt DFGLSbelt tconf tquantL tquantR Qconf;

zcrit = norminv(0.95,0,1);

%/* Choose lag length p by BIC */

p = bic(pred,maxp);

%/* Sample size */

n = size(ret,1);
T = n-p;

%/* Create variables */

r = ret(p+1:n);
x = pred(p+1:n);
Dx = pred(p+1:n)-pred(p:n-1);

Lx = [ones(T,1),pred(p:n-1)];

for i=1:p
    if i==1;
        Z = Lx;
    else
        Z = [Z,(pred(p-i+2:n-i+1)-pred(p-i+1:n-i))];
    end
end


%/* Estimate AR(1) for predictor variable */

coef = Lx\x;
v = x-Lx*coef;

Sv2 = v'*v/(T-2);

varmat = Sv2*inv(Lx'*Lx);
SErho = sqrt(varmat(2,2));

%/* Estimate AR(p) for predictor variable */

coef = Z\Dx;
e = Dx-Z*coef;

if p==1;
    psi = 0;
else
    psi = coef(3:p+1);
end

%/* Compute DF t-statistic */

varmat = (e'*e)/(T-p-1)*inv(Z'*Z);

df = coef(2)/sqrt(varmat(2,2));

%/* Compute DF-GLS statistic */

rhogls = 1-7/T;

xbar = [pred(p) ; pred(p+1:n)-rhogls*pred(p:n-1)];
inter = [1 ; (1-rhogls)*ones(T,1)];

mu = inter\xbar;

Zbar = Z(:,2)-mu;

if p>1
    Zbar = [Zbar,Z(:,(3:p+1))];
end

coef = Zbar\Dx; % why Dx not sth like Lx -mu?
resid = Dx-Zbar*coef;

varmat = (resid'*resid)/(T-p)*inv(Zbar'*Zbar);

dfgls = coef(1)/sqrt(varmat(1,1));

%/* Compute 95% confidence interval for c and rho */
conf = zeros(1,4);

%%% function Lconf in programs\theory\procs
conf(3) = Lconf(dfgls,DFGLSbelt,5);
conf(4) = Uconf(dfgls,DFGLSbelt,5);

conf(1) = 1+conf(3)/T;
conf(2) = 1+conf(4)/T;


%/* Predictive regression */

coef = Lx\r;
u = r-Lx*coef;

beta = coef(2);

%/* Estimate covariance matrix */

Sigma = ([u,e])'*([u,e])/(T-2);

delta = Sigma(1,2)/sqrt(Sigma(1,1)*Sigma(2,2));
norm = sqrt(Sigma(2,2)/Sigma(1,1));

omega2 = Sigma(2,2)/((1-sum(psi))^2);


%/* Estimate excess kurtosis */

Kurt = (mean([u,e].^4)./(mean([u,e].^2).^2))'-3*ones(2,1);
disp(Kurt');


%/* Compute point estimate and t-stat for beta */

normbeta = norm*beta;

varmat = Sigma(1,1)*inv(Lx'*Lx);
SEbeta = sqrt(varmat(2,2));

tstat = beta/SEbeta;

%/* Compute 90% confidence interval for beta using t-stat */

tbeta = zeros(1,2);


cL = BonLconf(df,DFbelt,-abs(delta),tconf(:,1));
cU = BonUconf(df,DFbelt,-abs(delta),tconf(:,1));
critL = BonLcrit(cL,cU,-abs(delta),tquantL);

cL = BonLconf(df,DFbelt,-abs(delta),tconf(:,2));
cU = BonUconf(df,DFbelt,-abs(delta),tconf(:,2));
critR = BonRcrit(cL,cU,-abs(delta),tquantR);

if delta<=0
    tbeta(1) = norm*(beta-critR*SEbeta);
    tbeta(2) = norm*(beta-critL*SEbeta);
else
    tbeta(1) = norm*(beta+critL*SEbeta);
    tbeta(2) = norm*(beta+critR*SEbeta);
end


% /* Compute 90% confidence interval for beta using Q-stat */

Qbeta = zeros(1,7);

cL = BonLconf(dfgls,DFGLSbelt,-abs(delta),Qconf(:,1));
cU = BonUconf(dfgls,DFGLSbelt,-abs(delta),Qconf(:,2));

rhoL = 1+cL/T;
rhoU = 1+cU/T;

correct = (T-2)/2*Sigma(1,2)/sqrt(Sigma(2,2)*omega2)*(omega2/Sv2-1)*SErho^2;

% /* At LOWER confidence interval for rho */

rbar = r-Sigma(1,2)/sqrt(Sigma(2,2)*omega2)*(x-rhoL*Lx(:,2));

coef = Lx\rbar;

Qbeta(1) = norm*(coef(2)+correct-zcrit*sqrt(1-delta^2)*SEbeta);
Qbeta(2) = norm*(coef(2)+correct+zcrit*sqrt(1-delta^2)*SEbeta);

% /* At UPPER confidence interval for rho */

rbar = r-Sigma(1,2)/sqrt(Sigma(2,2)*omega2)*(x-rhoU*Lx(:,2));

coef = Lx\rbar;

Qbeta(3) = norm*(coef(2)+correct-zcrit*sqrt(1-delta^2)*SEbeta);
Qbeta(4) = norm*(coef(2)+correct+zcrit*sqrt(1-delta^2)*SEbeta);

% /* Bonferroni confidence interval */

if delta<=0
    Qbeta(5) = Qbeta(3);
    Qbeta(6) = Qbeta(2);
else
    Qbeta(5) = Qbeta(1);
    Qbeta(6) = Qbeta(4);
end

% /* Restricted to rho<=1 */

if delta>0;
    Qbeta(7) = miss(1,1);
else
    rbar = r-Sigma(1,2)/sqrt(Sigma(2,2)*omega2)*Dx;

    coef = Lx\rbar;
    
    Qbeta(7) = norm*(coef(2)+correct-zcrit*sqrt(1-delta^2)*SEbeta);
end


return1=[p,delta,df,dfgls,conf];
return2=[tstat,normbeta,tbeta,Qbeta];


end
