/******************************************************
 * ML-estimation and LR-test of location for the para-*
 * meters of a  multivariate normal distribution with *
 * a monotone pattern of missing values.              *
 ******************************************************/

proc iml;

/*data from example 16.4.1*/
x={5 12  19,
   3  9 -21,
   5  7  14,
  -2  4  19,
   4 -2 -11,
   0  5  14,
  -6 15   4,
   1 -8   9,
   3  4   4,
   3  2 -31,
   0  1   .,
   3  7   .,
   4  8   .,
  -3 10   .,
   5 -1   .,
   5  .   .,
   5  .   .,
   3  .   .,
   5  .   .,
   5  .   .};

/*data from example 16.2.1 (inactive)*/
/*x={-8 -1,
    7  6,
   -2  4,
    0  2,
   -2  5,
    0  3,
   -2  .,
    1  .};*/
/*Sample size and number of variables:*/
n=nrow(x);
p=ncol(x);
print x n p;
xindex=((x^=.)#(1:p))[,<>]; /*vector of last observations for each
                               sample unit.*/
do i=1 to p;
ni=(xindex>=i)[+]; /*number of sample units that have an obs. on
                     the ith variable (and, because of monotony, on
                     every preceeding variable.)*/
Xi=X[1:ni,1:i];
xivec=Xi[,i];
if i=1 then do;
    muvec=Xi[:];
    Sigma=1/ni#(((xivec-muvec)##2)[+]);

    /*auxiliary calculations for the location test*/
    lambda=(Sigma/(1/ni#((xivec##2)[+])))##(ni/2);
    helprho=3/ni;
  end;
if i>1 then do;
   Xi=Xi[,1:(i-1)];
    /*auxiliary variables:*/
    comp=j(ni,1,1)||Xi;
    thetbet=inv(comp`*comp)*comp`*xivec; /*vector from (16.3.2))*/
    thetai=thetbet[1,];
    betai=thetbet[2:i,];
   /*calculation of mu:*/
     mui=thetai+betai`*muvec;
     muvec=muvec//mui;
   /*calculation of Sigma:*/
     sie=1/ni#xivec`*(I(ni)-comp*inv(comp`*comp)*comp`)*xivec;
     sigmaij=Sigma*betai;
     sigmaii=sie+betai`*Sigma*betai;
     Sigma=(Sigma||sigmaij)//(sigmaij`||sigmaii);
   /*calculations for the test of location:*/
     sietil=1/ni#xivec`*(I(ni)-Xi*inv(Xi`*Xi)*Xi`)*xivec;
     lambda=lambda#(sie/sietil)##(ni/2);
     helprho=helprho+(2#i+1)/ni;
  end;
end;

/*Test statistic for location test mu=0*/
  rho=1-helprho/(2#p);
  ll=-2#log(lambda);
  lladj=rho#ll;
  pvalue=1-probchi(lladj,p);
print "ML-estimates and location test for a single monotone sample";
print "Mean:" muvec;
print "Covariance:" Sigma;
print "Test of H0:mu=0:" lladj p pvalue;
quit;
run;