function [Psi,CV,CI,Nb,min_Nb,Bias,Skew,Var_psi]=matlab_pgm
%
%
% To execute the program, type
%
%    [Psi,CV,CI,Nb,min_Nb,Bias,Skew,Var_psi]=matlab_pgm
%
% at the Matlab prompt

load rock.dat
X=[ones(48,1) log(rock(:,2)/10000) log(rock(:,3)/2000) log(rock(:,4))];
Y=log(rock(:,5));
A=kron(eye(12),ones(4,1));
X=(1/4)*A'*X;
X(:,[2 3 4]')=exp(X(:,[2 3 4]'));
Y=(1/4)*A'*Y;
CC=[1 0.717 1.393 .203
    0 .320 .938 -0.132
    1 0 0 0
    0 1 0 0 
    0 0 1 0
    0 0 0 1];

ncc=size(CC,1);

[Psi,CV,CI,Nb,min_Nb,Bias,Skew,Var_psi]=Compute_ci(X,Y,CC');
Psi
CV{1}
CI{1}
Nb
min_Nb
Bias
Skew
Var_psi

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Psi,CV,CI,Nb,min_Nb,Bias,Skew,Var_psi]= ...
    Compute_ci(X,Y,C,conf1,conf2);
global c cx q0 A n
%
% Input Arguments:
% X: N x p model matrix.
%
% Y: N x 1 response vector.
%
% C: p x k matrix of coefficients.  The ith column of C contains
%    the p x 1 vector of coefficients for the ith linear function.
% conf1: confidence level for one-sided intervals.  Default is .95.
% conf2: confidence level for symmetric two-sided intervals.  Default is .95.
%
% Output arguments
% Psi: k x 1 vector of linear function estimates 
% CV: Cell array containing estimated critical values.
%     Specifically, CV{i} contains the estimated critical values
%     for the ith linear function.
% CI: Cell array containing the confidence interval endpoints.
%     Specifically, CV{i} contains the endpoints for the ith linear function.
% Nb: k x 1 vector of effective sample sizes.
% min_Nb: minimal effective sample size over all linear functions.
%
% Bias: k x 1 vector of estimates of rho_1(V)
%
% Skew: k x 1 vector of estimates of rho_3(V)
%
% Var_psi: : k x 1 vector of estimates Var(Psi)
%
% The CV{i} and CI{i} matrices contain critical values (estimated
%     percentiles of V) ordered according to Table 2 in the manuscript.

index=[1:12]';
Psi=[];
Var_psi=[];
Bias=[];
Skew=[];
Nb=[];
min_Nb=[];
if nargin < 3
  disp(['Error: Function Compute_ci requires a minimum of three'])
  disp(['input arguments.'])
  return
end
[p1,k]=size(C);
[N,p]=size(X);
if abs(p1-p)>0
  disp(['Error: The number of rows of C must match the number'])
  disp(['of columns of X'])
  return
end
n=N-rank(X);
adj=1/2;
XPX=X'*X;
XPXg=pinv(XPX);
Delta=C-XPX*XPXg*C;
Delta=reshape(Delta,p*k,1);
Delta=sqrt(Delta'*Delta);
if Delta > 1.e-10
  disp(['Error: Each column of C must lie in the row-space of X'])
  return
end
if nargin==3
  alpha1=.05;
  alpha2=.025;
elseif nargin==4
  alpha1=1-conf1;
  alpha2=.025;
else
  alpha1=1-conf1;
  alpha2=(1-conf2)/2;
end
alpha1=[alpha1 1-alpha1]';
alpha2=[alpha2 1-alpha2]';
crit0=tinv(alpha1,n);
crit_sym=tinv(alpha2,n);


A=speye(N)-X*XPXg*X';
min_Nb=1-min(diag(A));
min_Nb=1/min_Nb;
A3=A.^3;
A4=A.^4;
a=diag(A);
SSE=Y'*A*Y;
S=SSE/n;
Bhat=XPXg*X'*Y;
E=A*Y;
G2=sum(E.^2)/n;
G3=sum(E.^3)/n;
G4=sum(E.^4)/n;
G6=sum(E.^6)/n;
q2=a'*a/n;
q6=sum(sum(A3))/n;
q7=sum(sum(A4))/n;
dett=n/(q7*(n+2)-3*q2^2);
a4n=dett*(G4*(1+2/n-3*q2/n)+G2^2*3*(q7-q2));
a4=a4n/G2^2;
VarS2=max(0,dett*(q2*G4-(3*q2^2-2*q7)*G2^2)/n);
kappa3=(G3/(q6*G2^(3/2)))*(1+(3*VarS2)/(8*G2^2));
kappa4=max((a4n/G2^2)*(1+(VarS2)/(G2^2))-3,-2);
kappa32=min(kappa3^2,kappa4+2);

clear CV CI 
for ii=1:k
  CV{ii}=zeros(10,2);
  CI{ii}=zeros(10,2);
  c=C(:,ii);
  b=X*XPXg*c*sqrt(n);
  cx=b/sqrt(n);
  Nb_ii=b'*b/max(b.^2);
  Nb=[Nb; Nb_ii];
  q0=n*c'*XPXg*c;
  q1=a'*b/sqrt(n*q0);
  q3=sum(a.*(b.^2))/q0;
  q4=sqrt(n)*sum(b.^3)/(3*q0^(3/2));
  q5=n*sum(b.^4)/(3*q0^2);
  q8=sum(A3*b)/(sqrt(n*q0)*q6);
  q9=3*q5-q2-2*q3;
  q10=2*q3-q2-q5;
    psihat=c'*Bhat;
  Psi=[Psi; psihat];
  varhat=q0*S;
  Var_psi=[Var_psi; varhat];
  om1=-q1*kappa3/2;
  om12=q1^2*kappa32/4;
  om2=2+7*q1^2*kappa32/4+kappa4*(q2-q3);
  om2_star=2+kappa32*(-25*q1^2+36*q1*q4-10*q4^2)/4 ...
          +kappa4*(q2-q3+q8*(3*q1-2*q4));
  om3=3*kappa3*(q4-q1);
  Skew=[Skew; om3/sqrt(n)];
  Bias=[Bias; om1/sqrt(n)];
  om13=-3*q1*(q4-q1)*kappa32/2;
  om32=9*kappa32*(q4-q1)^2;
  om4=6+18*q1*(q1-q4)*kappa32+3*kappa4*(q5+q2-2*q3);
  om4_star=6-24*kappa32*(q1-q4)^2+3*kappa4*(q2+q5-2*q3+4*q8*(q1-q4));

  om2_star=n*(exp(om2_star/n)-1);
  om4_star=max(om4_star,-2*n);

  g7=kappa32*(q1^2+6*q1*q4-15*q4^2)/8+kappa4*q9/8;
  g8=kappa32*(q1-q4)*(q1-5*q4)/4+kappa4*q10/8;
  g9=om32/72+g8^2/(4*n);

  g10=(2+om4_star-4*om2_star)/8;
  g11=(6-om4_star)/24;
  g12=kappa32*q4*(q1-q4)/2-(2+om4_star-4*om2_star)/8;
  g13=(q4-q1)^2*kappa32/2-(6-om4_star)/24;
  g7n=max(exp(g7/n),adj/n^2);
  mg7n=max(exp(-g7/n),adj/n^2);
  g10n=max(exp(g10/n),adj/n^2);
  mg10n=max(exp(-g10/n),adj/n^2);
  g12n=max(exp(g12/n),adj/n^2);
  mg12n=max(exp(-g12/n),adj/n^2);

  d1=exp(-.5*(5-sqrt(17)))*((31-7*sqrt(17))*om3^2)/(72*n);
  
  m=zeros(6,1);
  m(2)=g7n;
  m(4)=g8/n;
  m(6)=-om32/(72*n);
  d2=solve_d(m,1/n);
  d3=0;
  if kappa3^2> 12*g11*g10n/(q1-q4)^2
    m=zeros(4,1);
    m(1)=kappa3*q4/(2*sqrt(n));
    m(2)=g10n;
    m(3)=kappa3*(q1-q4)/(2*sqrt(n));
    m(4)=g11/n;
    d3=solve_d(m,0.1);
  end
  d4=0;
  if kappa3^2> 12*g13*g12n/(q1-q4)^2
    m=zeros(4,1);
    m(1)=-kappa3*q4/(2*sqrt(n));
    m(2)=g12n;
    m(3)=-kappa3*(q1-q4)/(2*sqrt(n));
    m(4)=g13/n;
    d4=solve_d(m,0.1);
  end
  d5=max(0,-2*g11*exp(-.5)/(n*g10n));
  d6=max(0,2*g11*exp(-.5)/(n*mg10n));  

  CV{ii}(1,:)=crit0';
  tmp=(1+(om3/(2*sqrt(n)))*(om3/(6*sqrt(n))-om1/sqrt(n)-crit0));
  tmp=sign(tmp).*abs(tmp).^(1/3);
  Tcrit_Hall=(6*sqrt(n)/om3)*(1- tmp);
  CV{ii}(2,:)=Tcrit_Hall';
  Tcrit_expon2=solve_expon(crit0,crit0,d1, ...
              -om1/sqrt(n)+om3/(6*sqrt(n)),[1 -om3/(6*sqrt(n))]');
  CV{ii}(3,:)=Tcrit_expon2';
  Tcrit_expon3b=solve_expon(crit0,crit0,d3, ...
              kappa3*q4/(2*sqrt(n)), ...
             [g10n kappa3*(q1-q4)/(2*sqrt(n)) g11/n]');
  CV{ii}(4,:)=Tcrit_expon3b';
  Tcrit_3b=crit0.*g12n -q4*kappa3/(2*sqrt(n)) ...
        +exp(-d4*crit0.^2/2).*(crit0.^2*kappa3*(q4-q1)/(2*sqrt(n)) ...
                  +crit0.^3*g13/n);
  CV{ii}(5,:)=Tcrit_3b';
  Tcrit_expon3a=solve_expon(crit0,crit0,d5, ...
              0,[g10n 0 g11/n]');
  Tcrit_expon3=solve_expon(Tcrit_expon3a,Tcrit_expon3a,d1, ...
              -om1/sqrt(n)+om3/(6*sqrt(n)),[1 -om3/(6*sqrt(n))]');
  CV{ii}(6,:)=Tcrit_expon3';
  Tcrit_e3=crit0.*mg10n-exp(-d6*crit0.^2/2).*(crit0.^3*g11/n);
  Tcrit_3=solve_expon(Tcrit_e3,Tcrit_e3,d1, ...
              -om1/sqrt(n)+om3/(6*sqrt(n)),[1 -om3/(6*sqrt(n))]');
  CV{ii}(7,:)=Tcrit_3';
  CV{ii}(9,:)=crit_sym';
  Tcrit_expon_sym=solve_expon(crit_sym,crit_sym,d2,0, ...
             [g7n 0 g8/n 0 -om32/(72*n)]');
  CV{ii}(10,:)=Tcrit_expon_sym';
  Tcrit_sym=crit_sym.*exp(log(mg7n)-crit_sym.^2*g8/n+crit_sym.^4*g9/n);
  CV{ii}(11,:)=Tcrit_sym';

  rs=E./sqrt(a);
  rs=rs-mean(rs);
  Nboot=999;
%  Nboot=0;
%
% Set Nboot=0 to eliminate bootstrap intervals
%
  if Nboot > 0
    out_star0=boot_reg(rs);
    out_star=bootstrp(Nboot,'boot_reg',rs);
    out_star=[out_star; out_star0];
    Nboot=Nboot+1;
    s_hat=sqrt(S);
    boot_crit=prctile(out_star,100*alpha1);
    pct1=prctile(abs(out_star),100*alpha1(2));
    boot_crit_sym=[-sort(pct1); pct1];
    CV{ii}(8,:)=boot_crit';
    CV{ii}(12,:)=boot_crit_sym';
  else
    CV{ii}(8,:)=NaN(1,2);
    CV{ii}(12,:)=NaN(1,2);
  end
    
  CI{ii}=psihat-sqrt(varhat/n)*CV{ii};
  CI{ii}=CI{ii}(:,[2 1]');
  CV{ii}=[index CV{ii}];
  CI{ii}=[index CI{ii}];  
end

Nb=[[1:k]' Nb];
Psi=[[1:k]' Psi];
return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d=solve_d(mm,guess)
m0=mm(1);
k=length(mm)-1;
m=vec(mm(2:k+1));
d=guess;
del = 1;
while del > 1.e-8
  c=zeros(k+3);
  for j=2:k
    c(j-2+1)=c(j-2+1)+j*(j-1)*m(j);
    c(j+1)=c(j+1)-(2*j+1)*d*m(j);
    c(j+3)=c(j+3)+d^2*m(j);
  end
  c=c([k+3:-1:1]);
  r=roots(c);
  real_r=[];
  for j=1:length(r);
    if abs(imag(r(j))) < 1.e-10
      real_r=[real_r; real(r(j))];
    end
  end
  xd=real_r;
  c1=0;
  for j=2:k
    c1=c1+(j*xd.^(j-1)-d*xd.^(j+1))*m(j);
  end
  h1=m(1)+exp(-d*xd.^2/2).*c1;
  ii=find(h1==min(h1));
  xd=xd(ii(1));
  h1=h1(ii(1));
  s1=0;
  s2=0;
  for j=2:k
    s1=s1+((2*j+1)*xd^j-d*xd^(j+2))*m(j);
    s2=s2+(j*(j-1)*(j-2)*xd^(j-3)-d*(2*j+1)*xd^(j-1) ...
           +d^2*(j+2)*xd^(j+1))*m(j);
  end
  dxd=s1/s2;
  dh1=0;
  for j=2:k
    dh1=dh1+((j*(j-1)-d*(2*j+1)*xd^2+d^2*xd^4)*dxd-(j+2)*xd^3/2+d*xd^5/2) ...
       *xd^(j-2)*m(j);
  end
  del=max(-d,min(d/2,h1/dh1));
  d=d-del;
  del=abs(del);
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function BB=solve_expon(guess,AA,g,c0,c)

%
% Given vector AA, solve for vector BB, where
% A = c0 + B*c(1) + exp(-g*B^2/2)*sum_{i=2}^k(c_i*B^i)
%

B_out=[];
NA=length(AA);
c1=c(1);
c2=vec(c(2:end));
k=length(c);
pow=[2:k]';
for j=1:NA
  B=guess(j);
  A=AA(j);
  del=1;
  while del > 1.e-7
    h=c0-A+B*c1+exp(-g*B^2/2)*c2'*(B.^pow);
    H=c1+exp(-g*B^2/2)*(c2'*(pow.*(B.^(pow-1))) -g*B*c2'*(B.^pow));
    del=h/H;
    improve=-1;
    alpha=2;
    while improve<0
      alpha=alpha/2;
      d=del*alpha;
      Bnew=B-d;
      hnew=c0-A+Bnew*c1+exp(-g*Bnew^2/2)*c2'*(Bnew.^pow);
      improve=abs(h)-abs(hnew);
    end
    B=Bnew;
    del=abs(h);
  end
  B_out=[B_out; B];
end
BB=B_out;
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%