new;

      
   /*********************************************************
   **  Semiparametric estimation of a multiple-spell       **
   **  proportional hazards model                          **
   **********************************************************
   **  A Gauss program for Monte Carlo experiments         **
   **          [the censored case]                         **
   **********************************************************
   **      Estimation of the parametric part               ** 
   **      using Burke estimator                           **
   **********************************************************
   **      K_T: fourth order    K_X: second-order          **
   **********************************************************
   **      Censoring variable is observed                  ** 
   **********************************************************/



output file = beta_n1.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 ---@


   @--- Bandwidths ---@

if n==100;

        h1 = 4.5;
        h2 = h1;
        hg = 1.0;
        hx = 1.0;

elseif n==500;

        h1 = 3.5;
        h2 = h1;
        hg = 0.7;
        hx = 0.7;

endif;

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

       t2l = 3.5|0.5;
       x1l = 1|-1;
       x2l = 1|-1;


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

   @--- Bandwidths ---@

if n==100;

        h1 = 5.0;
        h2 = h1;
        hg = 1.2;
        hx = 1.2;

elseif n==500;

        h1 = 4;
        h2 = h1;
        hg = 0.9;
        hx = 0.9;

endif;

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

       t2l = 5|0.2;
       x1l = 1|-1;
       x2l = 1|-1;


endif;

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



    @--- Start of loop ---@ 


      seed = 12345; 
  result_b = {};
     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;
     alpha =  0.087;

   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] ---@ 


   /***************************************
   ***   Estimation of G(y1+y2|x1,x2)   ***
   ***************************************/ 

 @--- The kernel estimator of Pr(C > y1+y2|x1,x2) ---@


    y_c = y1 + y2; 
    pnx = kl2((x1 - x1')/hg).*kl2((x2 - x2')/hg); 
    pnx = meanc(pnx)/(hg^2);
    Hnc = (c .> y_c').*kl2((x1 - x1')/hg).*kl2((x2 - x2')/hg); 
   H_nc = meanc(Hnc)/(hg^2);
     Gn = H_nc./pnx;
     



   /***************************************
   ***    Estimation of beta            ***
   ***************************************/ 


    vmat = {};
      wx = (x1.<=x1l[1]).*(x1.>=x1l[2]).*(x2.<=x2l[1]).*(x2.>=x2l[2]);  @ trimming function @


       i = 1;
    do until i > n;

if wx[i]==0;    
     vhat = 1;
     vmat = vmat|vhat; 

else;

  x1e = x1[i];
  x2e = x2[i];
 
  @----- Quadrature method -----@      


    _intord = 3;
       vhat = intquad1(&ftn_b,t2l);   


     " rep - i - vhat - truev " rep~i~vhat~exp((x1e-x2e)*beta);  

@--- remove inaccurate estimates ---@

if (scalmiss(vhat)==1) or (abs(vhat)>1e+20) or (vhat<=0); 
 vhat = 1;
wx[i] = 0;
endif;

@-----------------------------------@
       vmat = vmat|vhat; 

endif;

     i = i + 1;
     endo;


    dt_x = wx.*(x1-x2);
     b_n = invpd(dt_x'(x1-x2))*(dt_x'ln(vmat));


result_b = result_b|b_n;

" estimator of beta " b_n; 


rep = rep + 1;
endo;


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


if (n==100) and (spec==1);       b1_100 = result_b;  save b1_100;
elseif (n==100) and (spec==2);   b2_100 = result_b;  save b2_100;
elseif (n==500) and (spec==1);   b1_500 = result_b;  save b1_500;
elseif (n==500) and (spec==2);   b2_500 = result_b;  save b2_500;
endif;


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



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


cls;

output on;
"============================================";
" the censored case ";
" specification " spec;
" sample size   " n;
" no. of rep.   "   n_rep;
" means of bandwidth h1 - h2 - hg - hx ";
h1~h2~hg~hx;
" rmse        " sqrt(meanc((result_b-beta)^2));
" mae         " median(abs(result_b-beta));
" mean bias   " meanc(result_b)-beta;
" std         " sqrt(meanc((result_b-meanc(result_b))^2));
" median bias " median(result_b)-beta;
output off;


   /********************************************************
   ***           Drawing the density of beta_n           ***
   ********************************************************/

dis = maxc(result_b)-minc(result_b);
bb = seqa(minc(result_b),dis/50,51);
hb = 1.06*stdc(result_b)*rows(result_b)^(-1/5);
fhat = kde(result_b,bb,hb);
library pgraph;
xy(bb,fhat);


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;



/**********************************************************/
/******    Integrand for the quadrature method    *********/
/**********************************************************/
proc ftn_b(t); 
local An,Bn,Rn,wt,f;

 @--- integrand ---@

    An = dt1.*dt2.*(y2 .> vecr(t)').*kl4((vecr(t)' - y1)/h1).*kl2((x1e - x1)/hx).*kl2((x2e - x2)/hx);
    Bn = dt1.*dt2.*(y1 .> vecr(t)').*kl4((vecr(t)' - y2)/h2).*kl2((x1e - x1)/hx).*kl2((x2e - x2)/hx);

    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])); 
     f = wt*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;

/**********************************************************/
/*************** Kernel Density Estimator *****************/
/**********************************************************/
proc(1) = kde(x,x_arg,h);

@
--- Inputs ---
    x = n*d vector of observations
x_arg = k*d vector of arguments of density
    h = d*1 vector of bandwidths
--- Output ---
k*1 vector of kernel density estimates
@

local n,d,k,i,arg,tmp,fhat;

   n = rows(x);
   d = cols(x);
   k = rows(x_arg);

 tmp = ones(k,n);
   i = 1;
   do until i > d;

          arg = (x_arg[.,i]-x[.,i]')/h[i];
          tmp = tmp.*(pdfn(arg)/h[i]);

   i = i + 1;
   endo;

fhat = meanc(tmp');

retp(fhat);
endp;