new;

/* Revised on April 2003 */
      
   /*********************************************************
   **  Semiparametric estimation of a multiple-spell       **
   **  proportional hazards model                          **
   **********************************************************
   **  A Gauss program for Monte Carlo experiments         **
   **          [the censored case]                         **
   **********************************************************
   **      Estimation of the nonparametric part            ** 
   **      using Burke estimator                           **
   **********************************************************/


output file = ld_n-c-r.out off;
format  /m1 /rd 12,4;
output off;

et = hsec;

  nset = {100, 500};
i_nset = 1;
do until i_nset > rows(nset);


         n = nset[i_nset];    @ sample size @
      beta = 1;               @ true beta @


   /****************************************************************
   ***         ln(Ld_0(t)) = -x*beta - u + ep                    ***
   *****************************************************************
   ***  Model specification 1: Weibull hazard                    *** 
   *****************************************************************
   ***  Baseline hazard function  : ld_0(t) = alpha*y            ***
   ***  Integrated baseline hazard: Ld_0(t) = alpha*0.5*y^2      ***
   *****************************************************************  
   ***  Model specification 2: U-shaped hazard                   *** 
   *****************************************************************
   ***   ld_0(t) = 0.05*(t/5)^(-2/3) + 0.57*(t/5)^5              ***
   ***   Ld_0(t) = 0.75*(t/5)^(1/3)  + 0.475*(t/5)^6             ***
   *****************************************************************/


specmat = {1,2};

i_spec = 1;
do until i_spec > rows(specmat);

spec = specmat[i_spec];


   /***************************************
   ***    Tuning Parameters             ***
   ***************************************/ 

      
if spec == 1; @--- model specification 1 ---@

if n==100;

        h1 = 3.5;
        h2 = h1;
        hz = 5.0;

elseif n==500;

        h1 = 2.5;
        h2 = h1;
        hz = 4.0;

endif;

    @--- Support of weighting functions ---@

       t2l = 3.5|0.5;
       z1l = 1|-1;
       z2l = 1|-1;

    @--- True ld_0 ---@

     alpha = 0.087;
      norm = ln(t2l[1]/t2l[2])/((t2l[1]-t2l[2])*alpha);  @ the integral of w(t)/ld(t) @

     t_min = 0.2;                                        @ the minimum of evaluation points @
     t_max = 3.8;                                        @ the maximum of evaluation points @
         m = 19;                                         @ the number of evaluation points @ 
     t1mat = seqa(t_min,0.2,m);                          @ evaluation points @
      ld_0 = alpha*t1mat;                                @ true ld_0 @ 

elseif spec == 2; @--- model specification 2 ---@

   @--- Bandwidths ---@

if n==100;

        h1 = 3.0;
        h2 = h1;
        hz = 6.0;

elseif n==500;

        h1 = 2.0;
        h2 = h1;
        hz = 4.0;

endif;
    @--- Support of the weighting functions ---@

       t2l = 5|0.2;
       z1l = 1|-1;
       z2l = 1|-1;

    @--- True ld_0 ---@

   _intord = 12;
      norm = intquad1(&ihazard,t2l[1]|t2l[2])/(t2l[1]-t2l[2]);   @ the integral of w(t)/ld(t) @ 

     t_min = 0.4;                                        @ the minimum of evaluation points @
     t_max = 5.5;                                        @ the maximum of evaluation points @
         m = 18;                                         @ the number of evaluation points @ 
     t1mat = seqa(t_min,0.3,m);                          @ evaluation points @
      ld_0 = 0.05*(t1mat/5)^(-2/3) + 0.57*(t1mat/5)^5;   @ true ld_0 @ 


endif;

   /***************************************
   ***    Monte Carlo experiments       ***
   ***************************************/ 



    @--- Start of loop ---@ 


      seed = 12345; 
    result = {};
     n_rep = 100;    @ the number of replications @
       rep = 1;            
   do until rep > n_rep;
   

   /***************************************
   ***    Generating data               ***
   ***************************************/ 

    x1 = rndns(n,1,seed);               @ covariates @
    x2 = rndn(n,1);
     u = (x1 + x2)/2;                   @ fixed effects @

    ep = ln(ln(1/(1-rndu(n,2))));       @ error term that has the CDF F(ep) =  1-exp(-exp(ep)). @

@====================================================================================@

if spec == 1;

   t_1 = sqrt((2/alpha)*exp(-x1*beta - u + ep[.,1]));  @ latent durations @
   t_2 = sqrt((2/alpha)*exp(-x2*beta - u + ep[.,2]));   

@====================================================================================@

elseif spec == 2;

    w1 = exp(-x1*beta - u + ep[.,1]);
    w2 = exp(-x2*beta - u + ep[.,2]);

   t_1 = {};
   t_2 = {};

     i = 1;
    do until i > n;

        coeff = zeros(19,1);
     coeff[1] = 0.475;
    coeff[18] = 0.75;
    coeff[19] = -w1[i]; 
          sol = polyroot(coeff); 
          sol = maxc(sol.*(real(sol).>0).*(imag(sol).==0)); @ pick up the correct root @
          t_1 = t_1|5*sol^3;     
      

    coeff[19] = -w2[i]; 
          sol = polyroot(coeff); 
          sol = maxc(sol.*(real(sol).>0).*(imag(sol).==0)); @ pick up the correct root @
          t_2 = t_2|5*sol^3;     
      

     i = i + 1;
     endo; 

@------------------------------------------------------------------ 
idea => To obtain the inverse of Ld_0, solve for x(=(t/5)^(1/3)):
              0.475*x^(18) + 0.75*x - w = 0.    
-------------------------------------------------------------------@


   t_1 = real(t_1);  @ observed t1 @
   t_2 = real(t_2);  @ observed t2 @
  
endif;

@====================================================================================@

if spec == 1;
   c_m = 1/20;              @ the inverse of the mean of censoring threshold @  

elseif spec == 2;
   c_m = 1/20;

endif;

     c = rndu(n,1);
     c = -(1/c_m)*ln(1-c); @ censoring threshold: exponential with mean 1/c_m @
    c1 = c;                 
   dt1 = (t_1 .< c1);
    c2 = (c1 - t_1).*dt1;
   dt2 = (t_2 .< c2); 

   y1 = minc((t_1~c1)');
   y2 = minc((t_2~c2)');

@--- Note: Approximately, meanc(dt1) = 0.78 and meanc(dt2) = 0.64 [spec. 1] 
                          meanc(dt1) = 0.87 and meanc(dt2) = 0.76 [spec. 2] ---@ 




   /***************************************
   ***  Loading the estimates of beta   ***
   ***************************************/ 

@--- Beta is estimated using the fourth-order kernel for K_T. --@

if (n==100) and (spec==1);
load b1_100;
  bhat = b1_100[rep];

elseif (n==100) and (spec==2);
load b2_100;
  bhat = b2_100[rep];

elseif (n==500) and (spec==1);
load b1_500;
  bhat = b1_500[rep];

elseif (n==500) and (spec==2);
load b2_500;
  bhat = b2_500[rep];

endif;


   z_1 = x1*bhat;
   z_2 = x2*bhat;

   /***************************************
   ***   Estimation of G(y1+y2)         ***
   ***************************************/ 

 @--- The estimator of Pr(C > y1+y2) ---@


    y_c = y1 + y2; 
    Hnc = (c .> y_c'); 
     Gn = meanc(Hnc);

     

   /********************************************************
   ***  Estimation of ld_0 by the quadrature method      ***
   ********************************************************/


   ldmat = {};

       i = 1;
    do until i > rows(t1mat);

        t1 = t1mat[i];


   @----- Quadrature method -----@      


    _intord = 3;
    _intrec = 0;
      ldhat = intquad3(&ftn_c,t2l,z1l,z2l)/norm; 

      ldmat = ldmat|ldhat; 

     " rep --- i " rep~i~ld_0[i]~ldhat; 

     i = i + 1;
     endo;

 result = result|ldmat';

rep = rep + 1;
endo;

@-------------------------------------------------------@


if (n==100) and (spec==1);       ld1_100 = result;  save ld1_100;
elseif (n==100) and (spec==2);   ld2_100 = result;  save ld2_100;
elseif (n==500) and (spec==1);   ld1_500 = result;  save ld1_500;
elseif (n==500) and (spec==2);   ld2_500 = result;  save ld2_500;
endif;


@-------------------------------------------------------@



library pgraph;
graphset;

 
    _pdate = "";
   _pltype = {3,6,6,6,6,6};

    xy(t1mat,ld_0~meanc(result));  
     xy(t1mat,ld_0~median(result)); 
     xy(t1mat,ld_0~(result[1:5,.]')); 



   /********************************************************
   ***           Printing results                        ***
   ********************************************************/

cls;
output on;
"============================================";
" the censored case ";
" specification " spec;
" sample size   " n;
" no. of rep.   "   n_rep;
" bandwidths h1 - h2 - hz ";
h1~h2~hz;
" integrated mse " meanc(meanc((result-ld_0')^2));
output off;

i_spec = i_spec + 1;
endo;

i_nset = i_nset  + 1;
endo;

et = hsec - et;
output on;
" Elapsed Time, in minutes "  et/6000;
output off;


end;



/**********************************************************/
/******    procedure to evaluate the integral     *********/
/**********************************************************/
proc (1) = ihazard(t);
local hf;
hf = 0.05.*((t./5).^(-2/3)) + 0.57.*((t./5).^5);  
hf = 1/hf;
retp(hf);
endp;


/**********************************************************/
/******    Integrand for the quadrature method    *********/
/**********************************************************/
proc ftn_c(t2,z1,z2); 
local An,Bn,Rn,wt,f;

@--- integrand ---@

    An = dt1.*dt2.*(y2 .> vecr(t2)').*kl2((vecr(t1)' - y1)/h1).*kl4((vecr(z1)' - z_1)/hz).*kl4((vecr(z2)' - z_2)/hz);
    Bn = dt1.*dt2.*(y1 .> vecr(t1)').*kl2((vecr(t2)' - y2)/h2).*kl4((vecr(z1)' - z_1)/hz).*kl4((vecr(z2)' - z_2)/hz);
    An = An./(Gn + 1e-10*(Gn.==0)); 
    Bn = Bn./(Gn + 1e-10*(Gn.==0)); 
 
    Rn = (sumc(An)/h1)./(sumc(Bn)/h2);  
    wt = (1/(t2l[1]-t2l[2]))*kl2(vecr(z1)).*kl2(vecr(z2)); 
     f = wt.*exp(vecr(z2)-vecr(z1)).*Rn;  
     f = f';  

retp(f);
endp;



/**********************************************************/
/*************** 2nd Order Kernel Function ****************/
/**********************************************************/
proc(1) = kl2(v);
local term;
term = (15/16)*((1 - v^2)^2);
term = term.*(abs(v) .<= 1);
retp(term);
endp;


/**********************************************************/
/*************** 4th Order Kernel Function ****************/
/**********************************************************/
proc(1) = kl4(v);
local term;
term = (105/64)*(1 - 5*v^2 + 7*v^4 - 3*v^6);
term = term.*(abs(v) .<= 1);
retp(term);
endp;


/**********************************************************/
/*************** 6th Order Kernel Function ****************/
/**********************************************************/
proc(1) = kl6(v);
local term;
term = (315/2048)*(15 - 140*v^2 + 378*v^4 - 396*v^6 + 143*v^8);
term = term.*(abs(v) .<= 1);
retp(term);
endp;