// ---------------------------------------------------------------------
//      				Solution and estimation of a 
//         				  Model of Cow Replacement
//
//
//
//                  Victor Aguirregabiria and Arvind Magesan
//                             April 8, 2013
//
// ---------------------------------------------------------------------

new;
cls;
library pgraph;
seed =2393457892;

//
// -----------------------------------------------------------------
// MILOGIT - Estimation of a Logit Model by Maximum Likelihood
//           The optimization algorithm is a Newton's method.
//
// by Victor Aguirregabiria
//
// Last revision: September 2001
// -----------------------------------------------------------------
//
// Format      {best,varest} = milogit(ydum,x,namesb,inibeta,tol)
//
// Input        
//	ydum	- vector of observations of the dependent variable
//    x	- matrix of explanatory variables
//	zrest	- vector of explanatory variable restricted to have beta = 1
//	namesb- vector with names of parameters
//	tol	- scalar for the level of tolerance of convergence 
//
// Output       best    - ML estimates
//              varest  - estimate of the covariance matrix
// -----------------------------------------------------------------
//

proc (1) = loglogit(ydum,x,zr,b) ;
  local expxb, Fxb, llik, myzero ;
  myzero = 1E-12 ;
  expxb = exp(-x*b-zr) ;
  Fxb = 1./(1+expxb) ;
  Fxb = Fxb + (myzero - Fxb).*(Fxb.<myzero)
            + (1-myzero - Fxb).*(Fxb.>1-myzero);
  llik = ydum'*ln(Fxb) + (1-ydum)'*ln(1-Fxb) ;
  retp(llik) ;
endp ;


proc (2) = milogit(ydum,x,zrest,namesb,inibeta,tol) ;
  local nobs, nparam, eps1, eps2, iter, criter1, 
        criter2, expxb0, Fxb0, phixb0, lamdab0, 
        dlogLb0, d2logLb0, b0, b1, lamda0, lamda1, 
        Avarb, sdb, tstat, llike, numy1, numy0, 
        logL0, LRI, pseudoR2, k ;

  format /mb1 /ros 16,6 ;

  nobs = rows(ydum) ;
  nparam = cols(x) ;
  eps1 = tol ;
  eps2 = tol ;
  b0 = inibeta ;  
  iter=1 ;
  criter1 = 1000 ;
  criter2 = 1000 ;

  do while (criter1>eps1).OR(criter2>eps2) ;
/*
    "" ;
    "Iteration                = " iter ;
    "Log-Likelihood function  = " loglogit(ydum,x,b0) ;
    "Norm of b(k)-b(k-1)      = " criter1 ;
    "Norm of Gradient         = " criter2 ;
    "" ;
*/    
    expxb0 = exp(-x*b0-zrest) ;
    Fxb0 = 1./(1+expxb0) ;
    dlogLb0 = x'*(ydum - Fxb0) ;
    d2logLb0 =  ( (Fxb0.*(1-Fxb0)).*x )'*x ;
    b1 = b0 + invpd(d2logLb0)*dlogLb0 ;
    criter1 = sqrt( (b1-b0)'*(b1-b0) ) ;
    criter2 = sqrt( dlogLb0'dlogLb0 ) ;
    b0 = b1 ;
    iter = iter + 1 ;
  endo ;

  expxb0 = exp(-x*b0-zrest) ;
  Fxb0 = 1./(1+expxb0) ;
  Avarb = - ( (Fxb0.*(1-Fxb0)).*x )'*x ;
  Avarb = inv(-Avarb) ;
  sdb    = sqrt(diag(Avarb)) ;
  tstat  = b0./sdb ;
  llike  = loglogit(ydum,x,zrest,b0) ;
  numy1  = sumc(ydum) ;
  numy0  = nobs - numy1 ;
  logL0  = numy1*ln(numy1) + numy0*ln(numy0) - nobs*ln(nobs) ;
  LRI    = 1 - llike/logL0 ;
  pseudoR2 = 1 - ( (ydum - Fxb0)'*(ydum - Fxb0) )/numy1 ;

  "Number of Iterations     = " iter ;
  "Log-Likelihood function  = " llike ;
  "Likelihood Ratio Index   = " LRI ;
  "Pseudo-R2                = " pseudoR2 ;
  "" ;
 "     ----------------------------------------------------------------";
 "         Parameter       Estimate        Standard        t-ratios";
 "                                         Errors" ;
 "     ----------------------------------------------------------------";
  k=1;
  do while k<=nparam;
    print $namesb[k];;b0[k];;sdb[k];;tstat[k];
    k=k+1 ;
  endo;
 "     ----------------------------------------------------------------";

  retp(b0,avarb) ;
endp ;


//	---------------------
//		0. Constants
//	---------------------
beta = 0.95 ;   // Discount factor
nomega = 201 ;  // Number of cells in discretization of omega
maxlife = 5 ;   // Maximum number of lactations
numsdome = 5 ;  // Number of standard deviations for minimum and maximum values in the grid of omega
agestate = seqa(1,1,maxlife); // Grid of possible cow ages - possible values of y_it


//	-----------------------
//		1. Reading data
//	-----------------------
nobs = 2340 ;
nvar = 6 ;
load x[nobs,nvar] = 
c:\MYPAPERS\GENERALIZED_EULER_EQUATIONS\data\mir-sch.dat ;
farm   = x[.,1] ;			//Farm ID
profit = x[.,2] ;			//Profit
milk   = x[.,3] ;			//Milk 
rcost  = x[.,4] ;			//Replacement Cost ct
lact   = x[.,5] ;			//Age
drep   = x[.,6] ;			//Replacement Decision
clear x ;


// Creates cow identifier
cow = zeros(nobs,1) ;
cow[1] = 1 ;
i=2 ;
do while i<=nobs ;
  cow[i] = cow[i-1] + (lact[i]==1) ;
  i=i+1 ;
endo ;numcows = maxc(cow);



//	----------------------------------------
//		2. Selection of complete histories
//	----------------------------------------
//	this keeps cows for whom we observe their entire history

select = zeros(nobs,1) ;
j=1 ;
ind2 = 0 ;
do while j<=numcows ;
  buff = selif(drep,cow.==j) ;
  numj = rows(buff) ;
  ind1 = ind2 + 1 ;
  ind2 = ind1 + numj - 1 ;
  select[ind1:ind2] = sumc(buff)*ones(numj,1) ;
  j=j+1 ;
endo ;

farm   = selif(farm,select) ;
profit = selif(profit,select) ;
milk   = selif(milk,select) ;
rcost  = selif(rcost,select) ;
lact   = selif(lact,select) ;
drep   = selif(drep,select) ;
cow    = selif(cow,select) ;

clear select ;

nobs = rows(milk) ;

//Creates lag of milk production
milk_1 = (lact.>1).*(0 | milk[1:nobs-1]) ;

//	----------------------------------------
//		3. Descriptive Statistics    
//	----------------------------------------
//	a) The probability of replacement as a function of age
prep1 = ((lact.==1)'drep)/sumc((lact.==1));
prep2 = ((lact.==2)'drep)/sumc((lact.==2));
prep3 = ((lact.==3)'drep)/sumc((lact.==3));
prep4 = ((lact.==4)'drep)/sumc((lact.==4));
prep5 = ((lact.==5)'drep)/sumc((lact.==5));

prep = meanc(drep);
vrep = vcx(drep);
preptab = prep1|prep2|prep3|prep4|prep5|prep|vrep;

//	b) Milk production as a function of age
emilk1 = ((lact.==1)'milk)/sumc((lact.==1));
emilk2 = ((lact.==2)'milk)/sumc((lact.==2));
emilk3 = ((lact.==3)'milk)/sumc((lact.==3));
emilk4 = ((lact.==4)'milk)/sumc((lact.==4));
emilk5 = ((lact.==5)'milk)/sumc((lact.==5));

emilk = meanc(milk);
vmilk = vcx(milk);
emilktab = emilk1|emilk2|emilk3|emilk4|emilk5|emilk|vmilk;
histp(selif(lact,drep.==1),100);

//	----------------------------------------
//		4. Estimation of production function
//	----------------------------------------
//	Production function has the form:
//		ln(m_{it}) = alpha_1 * 1{y_{it}=1} + ... 
//			+ alpha_5 * 1{y_{it}=5} + omega_{it}
//
//	where: omega_{it} = rho * omega_{it-1} + psi_{it}
//	----------------------------------------

vage = seqa(1,1,maxlife) ;
dumlact = (lact.==1)~(lact.==2)~(lact.==3)~(lact.==4)~(lact.==5) ;
lnmilk = ln(milk);
lnmilk_1 = (lact.>1).*ln((lact.==1) + milk_1); 

// lnmilk_1 is set at 0 for the lact==1								// This is helpful to construct the FGLS

//	Estimation proceeds in 3 steps. 

//	Step 1: 
//	First we estimate alpha_1 using the fact 
//	that alpha_1 = E[ln(m_{it})| y_{it}=1]

alp1hat = sumc(lnmilk.*(lact.==1))/sumc(lact.==1);

//	Can export data to STATA to estimate standard errors etc.
statadata = lnmilk~dumlact;

//	redsamp1 = indexcat(dumlact[.,1], maxc(dumlact[.,1]));
//	yvar = lnmilk;
//	xvar = dumlact[.,1] ;
//	__con=0;
//	{vnam,m,best,stb,vc,stderr,sigma,cx,rsq,resid,dwstat} = ols(0,yvar,xvar) ;


//	Step 2: 
//	We use the structure of omega_{it} to estimate:
//		ln(m_{it}) = rho * ln(m_{it-1}) + gamma_2 * 1{y_{it}=2} +
//			 ... + gamma_5 * 1{y_{it}=5} + psi_{it}
//	where gamma_j = alpha_j - rho * alpha_{j-1}

redsampi = indexcat(dumlact[.,1], minc(dumlact[.,1]));
redsamp = lnmilk[redsampi]~lnmilk_1[redsampi]~dumlact[redsampi,2:cols(dumlact)] ;

"" ;
"--------------------------------------------------------" ;
"Production function estimation controlling for selection" ;
"--------------------------------------------------------" ;
"" ;
explan = redsamp[.,2:cols(redsamp)];
yvar = redsamp[.,1];
__con = 0 ;
__altnam = "milk_1"|"lact2"|"lact3"|"lact4"|"lact5"|"milk" ;
{vnam,m,best,stb,vc,stderr,sigma,cx,rsq,resid,dwstat} = ols(0,yvar,explan) ;
rhohat = best[1];

//	Step 3: Give the estimate of rho, we can get the FGLS estimator
//			of alpha1, alpha2, ..., alpha5.
yvar = lnmilk - rhohat * lnmilk_1 ;
explan = ((lact.==1) - rhohat.*(lact.==2))	   
       	~ ((lact.==2) - rhohat.*(lact.==3))
	~ ((lact.==3) - rhohat.*(lact.==4))	   
	~ ((lact.==4) - rhohat.*(lact.==5))	   
	~ ((lact.==5));
__con = 0 ;

__altnam = "lact1*"|"lact2*"|"lact3*"|"lact4*"|"lact5*"|"milk*" ;
{vnam,m,best,stb,vc,stderr,sigma,cx,rsq,resid,dwstat} = ols(0,yvar,explan) ;

alphat = best;
omega = lnmilk - ((lact.==1)~(lact.==2)~(lact.==3)~(lact.==4)~(lact.==5))*alphat;
	   
//	To obtain the transition of omega we use:
//	omega_{it} = rho * omega_{it-1} + xi_{it}

xi = omega-rhohat*(lact.>1).*(0 | omega[1:nobs-1]) ;
sigxi = stdc(xi);

//	----------------------------------------
//		5. Discretization omega
//	----------------------------------------
sdomega = stdc(omega);
minome = -numsdome * sdomega ;
maxome = numsdome * sdomega ;
stepome = int(10000*(maxome-minome)/(nomega-1))/10000 ;
omeval = seqa(minome,stepome,nomega) ;

//	----------------------------------------
//		6. Transition probabilities of omega
//	----------------------------------------
//	A column of ometran contains the whole probability distribution
//	of omega_{t+1} conditional on a particular value of omega_{t}

ometh0 = -1e12 | (omeval[1:nomega-1]+stepome) ;
ometh1 = (omeval[1:nomega-1]+stepome) | 1e12 ;
ometran = cdfn((ometh1-rhohat*omeval')/sigxi) - cdfn((ometh0-rhohat*omeval')/sigxi) ;
ometran = ometran./(sumc(ometran)') ;
ometran0 = ometran ;
ometran1 = ometran[.,round(nomega/2)].*.ones(1,nomega) ;
clear ometran ;

//Discretize true values of omega to the grid 

omfind = abs(omega- (omeval').*ones(rows(omega),1)  ) ;
omfind = omfind';
omfound = {};
i0 = 1;
do while i0 <= cols(omfind);
  omfound = omfound|indexcat(omfind[.,i0], minc(omfind[.,i0]));
  i0 = i0+1;
endo;

omega_d = omeval[omfound];xy(omega,omega~omega_d) ;

//	----------------------------------------
//		7. Transition age
//	----------------------------------------
//	A column of agetran0 contains the whole probability 
//	distribution (degenerate)
//	of age_{t+1} conditional on a particular value 
//	of age_{t} and no replacement at t.
//	Similarly, a column of agetran1 contains 
//	the whole probability distribution (degenerate)
//	of age_{t+1} conditional on a particular value 
//	of age_{t} and replacement at t.

agetran0 = seqa(2,1,maxlife-1) | 1 ;
agetrans0 = zeros(maxlife,maxlife);
i0 = 1;
do while i0<=maxlife;
  agetrans0[agetran0[i0],i0] = 1;
  i0 = i0+1;
endo;

agetran1 = ones(maxlife,1) ;
agetrans1 = zeros(maxlife,maxlife);
i0 = 1;
do while i0<=maxlife;
agetrans1[agetran1[i0],i0] = 1;
i0 = i0+1;
endo;

//	---------------------------------------------
//		8. Transition Price (profit contribution)
//	---------------------------------------------
//	One option is to assume it is constant and equal 
//	to the mean of the price 

pt = meanc(profit);


//	---------------------------------------------
//		9. Replacement cost Transition    
//	---------------------------------------------
//	One option is to assume it is constant and equal 
//	to the mean of the cost

ct = meanc(rcost);
 
//	---------------------------------------------
//		10. Estimation                     
//	---------------------------------------------
states = (agestate.*.ones(nomega,1))~(ones(maxlife,1).*.omeval);
Nstates = rows(states);
agepay = states[.,1].==(agestate');

//	---------------------------------------------
//		10.(a) TWO-STEP PML 
//	---------------------------------------------
//	Step 1.	Estimation of replacement probabilities 
//	by Victor's Logit code

agedat = (lact.==1)~(lact.==2)~(lact.==3)~(lact.==4);
dat = agedat~(agedat.*omega_d)~(agedat.*(omega_d.^2));

//Remove age=5 observations
tooold = (lact.==5);
tooold = indexcat(tooold, minc(tooold));

dat = dat[tooold,.];
drep = drep[tooold];

pnames = "age1"|"age2"|"age3"|"age4"       
	|"age1om"|"age2om"|"age3om"
	|"age4om"|"age1omsq"|"age2omsq"
	|"age3omsq"|"age4omsq";

tolerance = 1e-4;
inibeta = zeros(rows(pnames),1);
restz = zeros(rows(drep),1);

{brd_logit,varest} = milogit(drep,dat,restz,pnames,inibeta,tolerance);

//	Using reduced form estimates to construct 
//	initial CPs

redstat = agepay[.,1:4]
	~(agepay[.,1:4].*states[.,2])
	~(agepay[.,1:4].*(states[.,2].^2));
initcp = redstat*brd_logit;
initcp= 1/(1+exp(-initcp));

//Picture of replacement probability by age
age1 = initcp[1:nomega];
age2 = initcp[nomega+1:2*nomega];
age3 = initcp[2*nomega+1:3*nomega];
age4 = initcp[3*nomega+1:4*nomega];
graphset ;
_plwidth = 20 ;
_pdate = "" ;
xlabel("Omega") ;
ylabel("Replacement Probability") ;
_plegstr = "1 Year\0002 Year\0003 Year\0004 Year" ;
_plegctl = (2|4|2|2) ;
xy(omeval,age1~age2~age3~age4);

// Assign a state ID to each data point
state_id = {};
d0 = 1;
real_data = lact~omega_d;
do while d0<=nobs;
	sfind = sumr(real_data[d0,.].==states).==cols(states);
	sfind = indexcat(sfind,maxc(sfind));
	state_id = state_id|sfind;
	d0 = d0+1;
endo;

// Reduced state_ids not including age 5
tooold = states[.,1].==5;
tooold = indexcat(tooold, maxc(tooold));
tooold = minc(tooold);
young = 1-(state_id.>=tooold);
young = indexcat(young, maxc(young));

state_idred = state_id[young];

//Milk production with discretized omega 
MilkRev_d = exp(agepay*alphat + states[.,2]);

//Transition Matrices 
F0 = agetrans0.*.ometran0;
F1 = agetrans1.*.ometran1;
F0 = F0';
F1 = F1';

CPe_0 = initcp;
CPe_0  = (states[.,1].<=4).*CPe_0  
	+ (states[.,1].==5)*1; 

//	Setting CCPS at age 5 to 1
//	Payoffs with age dummies in cost
//	pay0 = (pt*MilkRev_d)~(-agepay[.,3:cols(agepay)])~(zeros(rows(agepay),1));
//	pay1 = (pt*MilkRev_d)~(-agepay[.,3:cols(agepay)])~(-ones(rows(agepay),1));

//Payoffs with linear functions of age

pay0 = (pt*MilkRev_d)
	~(-states[.,1])~zeros(nstates,1);
pay1 = (pt*MilkRev_d)
	~(-states[.,1])
	~(-ones(nstates,1));
z0 = (states[.,1].<=4).*(pay0)
   + (states[.,1].==5).*(pay1);
z1 = pay1;
e0 = 0.5772156649 - ln((states[.,1].==5) 
	+ (states[.,1].<=4).*(1-CPe_0));
e1 = 0.5772156649 - ln((states[.,1].==5) 
	+ (states[.,1].<=4).*CPe_0);
zbar_P = CPe_0.*z1 + (1-CPe_0).*z0;
ebar_P = (CPe_0.*e1 + (1-CPe_0).*e0);

Fbar_P = CPe_0.*F1 + (1-CPe_0).*F0;
W_Pz = inv(eye(nstates) - beta*Fbar_P)*zbar_P;
W_Pe = inv(eye(nstates) - beta*Fbar_P)*ebar_P;
ztilde0 = z0 + beta*F0*W_Pz;
ztilde1 = z1 + beta*F1*W_Pz;
etilde0 =  beta*F0*W_Pe;
etilde1 =  beta*F1*W_Pe;

ztilde = (ztilde1[state_idred,.]-ztilde0[state_idred,.]);
etilde = etilde1[state_idred,.]-etilde0[state_idred,.];
Phat1 = CPe_0[state_idred]; 
oddsrat = ln(Phat1./(1-Phat1));
twosteppml = inv(ztilde'ztilde)*(ztilde'(oddsrat-etilde));

STatadata = (oddsrat-etilde)~ztilde;
yvar = ln(Phat1./(1-Phat1))-etilde; 
explan = ztilde;
__altnam = "1/sigma"|"MC/sigma"|"RC/sigma"|"depvar" ;
{vnam,m,best,stb,vc,stderr,sigma,cx,rsq,resid,dwstat} = ols(0,yvar,explan) ;
pnames = "1/sigma"|"MC/sigma"|"RC/sigma";       
tolerance = 1e-4;
yvar = drep; 
explan = ztilde;
restz = etilde;
inibeta = best;
{pml_1,var_pml_1} = milogit(drep,explan,restz,pnames,inibeta,tolerance);
sigmahat = 1/pml_1[1];
mchat = pml_1[2]/pml_1[1];
rchat = pml_1[3]/pml_1[1];
jacobian =  ((-sigmahat) ~ 0 ~ 0)		 
	|  ((-mchat)    ~ 1 ~ 0)		 
	|  ((-rchat)    ~ 0 ~ 1);
jacobian =  sigmahat * jacobian;
vctheta = jacobian'*(var_pml_1*jacobian);
setheta = sqrt(diag(vctheta));  

"" ;  
"Number of Observations   = " rows(yvar);
"Pseudo-R2                = " rsq ;
"" ;
"     ----------------------------------------------------------------";
"         Parameter       Estimate        Standard        t-ratios";
"                                         Errors" ;
"     ----------------------------------------------------------------";  
"	SIGMA	";;	sigmahat;; setheta[1];; (sigmahat/setheta[1]);  
"	MC   	";;	mchat;;    setheta[2];; (mchat/setheta[2]);
"	RC   	";;	rchat;;    setheta[3];; (rchat/setheta[3]);
"     ----------------------------------------------------------------";

//	---------------------------------------------
//		10.(b) NPL 
//	---------------------------------------------
//	Now iterating until convergence in the 
//	NPL mapping - yields the MLE.

tol = 1e-10;
CPe_0 = initcp;
CPe_1 = ones(nstates,1);
 
crit = maxc(abs(CPe_1-CPe_0)) ;
it = 0;
do while (crit>=tol);  
  "NPL ITERATION	=";; it;
  CPe_0  = (states[.,1].<=4).*CPe_0  + (states[.,1].==5)*1;

  e0 = 0.5772156649 - ln((states[.,1].==5) 
	+ (states[.,1].<=4).*(1-CPe_0));  
  e1 = 0.5772156649 - ln((states[.,1].==5) 
	+ (states[.,1].<=4).*CPe_0);

  zbar_P = CPe_0.*z1 + (1-CPe_0).*z0;
  ebar_P = (CPe_0.*e1 + (1-CPe_0).*e0);

  Fbar_P = CPe_0.*F1 + (1-CPe_0).*F0;
  W_Pz = inv(eye(nstates) - beta*Fbar_P)*zbar_P;
  W_Pe = inv(eye(nstates) - beta*Fbar_P)*ebar_P;
  ztilde0 = z0 + beta*F0*W_Pz;
  ztilde1 = z1 + beta*F1*W_Pz;
  etilde0 =  beta*F0*W_Pe;
  etilde1 =  beta*F1*W_Pe;

  ztilde = (ztilde1[state_idred,.]-ztilde0[state_idred,.]);
  etilde = etilde1[state_idred,.]-etilde0[state_idred,.];
  Phat1 = CPe_0[state_idred]; 
  oddsrat = ln(Phat1./(1-Phat1));
  NPLest = inv(ztilde'ztilde)*(ztilde'(oddsrat-etilde));
  pnames = "1/sigma"|"MC/sigma"|"RC/sigma";         
  tolerance = 1e-4;  
  yvar = drep;   
  explan = ztilde;  
  restz = etilde;  
  inibeta = NPLest;
  {NPLest,var_NPLest} = milogit(drep,explan,restz,pnames,inibeta,tolerance);  
  ztilde0u = (z0+ beta*F0*W_Pz)*NPLest ;
  ztilde1u = (z1+ beta*F1*W_Pz)*NPLest ;

  kern = (ztilde1u-ztilde0u) + (etilde1-etilde0) ;
  CPe_1= 1/(1+exp(-kern));  

  //  Forcing CCPS at age 5 to be 1
  CPe_1 = (states[.,1].<=4).*CPe_1  + (states[.,1].==5)*1;
  crit = maxc(abs(CPe_1-CPe_0)) ;
  CPe_0 = CPe_1;
  it = it+1;
endo;

sigmahat = 1/NPLest[1];
mchat = NPLest[2]/NPLest[1];
rchat = NPLest[3]/NPLest[1];
jacobian =  ((-sigmahat) ~ 0 ~ 0)		 
	|  ((-mchat)    ~ 1 ~ 0)			|  ((-rchat)    ~ 0 ~ 1);
jacobian =  sigmahat * jacobian;		 
vctheta = jacobian'*(var_NPLest*jacobian);
setheta = sqrt(diag(vctheta));  

"" ;  
"Number of Observations   = " rows(yvar);
"Pseudo-R2                = " rsq ;
"" ;
"     ----------------------------------------------------------------";
"         Parameter       Estimate        Standard        t-ratios";
"                                         Errors" ;
"     ----------------------------------------------------------------";  
"	SIGMA	";;	sigmahat;; setheta[1];; (sigmahat/setheta[1]);  
"	MC   	";;	mchat;;    setheta[2];; (mchat/setheta[2]);
"	RC   	";;	rchat;;    setheta[3];; (rchat/setheta[3]);
"     ----------------------------------------------------------------"; 

//Picture of replacement probability by age
age1 = CPe_0[1:nomega];
age2 = CPe_0[nomega+1:2*nomega];
age3 = CPe_0[2*nomega+1:3*nomega];
age4 = CPe_0[3*nomega+1:4*nomega];
graphset ;
_plwidth = 20 ;
_pdate = "" ;
xlabel("Omega") ;
ylabel("Replacement Probability") ;
_plegstr = "1 Year\0002 Year\0003 Year\0004 Year" ;
_plegctl = (2|4|2|2) ;
_pmcolor = {0, 0, 0, 0, 0, 0, 0, 0, 15 };
xy(omeval,age1~age2~age3~age4);

//	-------------------------------------------------------
//		10.(c)	Estimation by Euler Equations - GMM
//	-------------------------------------------------------
//  The Euler equation is:
//  E_{t}[ M_{t+1}- theta_{R}+ theta_{C}*y_{t}+ sigma*e_{t+1} ]=0
//
//	where
//	M_{t+1} = beta * [M(1,omega_{t+1})- M(y_{t}+1,omega_{t+1})]
//		e_{t+1} = [lnP(0|x_{t})+ lnP(1|y_{t}+1,omega_{t+1})]
//	- [lnP(1|x_{t})+ lnP(1|1,omega_{t+1}) ]
//	-------------------------------------------------------

// Random draws of psi_{t+1}
drawxi_0 = int(nobs*rndu(nobs,1));
drawxi_0 = drawxi_0 + (drawxi_0.==0);
drawxi_0 = xi[drawxi_0];
drawxi_1 = int(nobs*rndu(nobs,1));
drawxi_1 = drawxi_1 + (drawxi_1.==0);
drawxi_1 = xi[drawxi_1];

// Construction of M_{t+1} (mtilde)
omegatod =  omega;    

// Today's omega
omegatom_0 = rhohat*omegatod + drawxi_0;	

// Tomorrow's omega if NO replacement today
omegatom_1 = drawxi_1; 	

// Tomorrow's omega if replacement today
ytod = lact;          	

// Today's age
ytodp1 = ytod+1;        

// Today's age plus 1
alpha_yp1 = (ytodp1.<=5).*ytodp1 + (ytodp1.>=6).*1;
alpha_yp1 = alphat[alpha_yp1];

mtilde = beta*(exp(alphat[1] + omegatom_0) 
	- exp(alpha_yp1 + omegatom_0));

// Construction of M_{t} (mdiff)
mdiff = beta*(exp(alphat[1] + omegatod) 
	- exp(alphat[ytod] + omegatod));

// Construction of e_{t+1} (etilde)
//		Today's ccps 
lackern = (ytod.==1)~(ytod.==2)~(ytod.==3)~(ytod.==4);
kern = lackern~(lackern.*omegatod)~(lackern.*(omegatod.^2));
kern = kern*brd_logit;
phatp = 1/(1+exp(-kern));
phatp = (ytod.<=4).*phatp + (ytod.==5);

//	Tomorrow's ccps p(1, omega_{t+1})
lackern = ones(nobs,1)~zeros(nobs,3);
kern = lackern~(lackern.*omegatom_1)~(lackern.*(omegatom_1.^2));
kern = kern*brd_logit;phatp1_1 = 1/(1+exp(-kern));

//	Tomorrow's ccps p(y_{t}+1, omega_{t+1})
lackern = (ytodp1.==1)~(ytodp1.==2)~(ytodp1.==3)~(ytodp1.==4);
kern = lackern~(lackern.*omegatom_0)~(lackern.*(omegatom_0.^2));
kern = kern*brd_logit;
phatp_1y = 1/(1+exp(-kern));
phatp_1y = (ytodp1.<=4).*phatp_1y + (ytodp1.==5).*1;
etilde = (ln(1-phatp) + beta*ln(phatp_1y)) 	   
	- (ln(phatp)   + beta*ln(phatp1_1));

// IV Regressionyobs = selif(mtilde,ytod.<=4);
xobs = ones(nobs,1)~(-beta*ytod)~(-beta*etilde);
xobs = selif(xobs,ytod.<=4);
zinst = ones(nobs,1)~ytod~omegatod~ln(phatp)~mdiff;
zinst = selif(zinst,ytod.<=4);
xbuff = inv(zinst'*zinst)*(zinst'*xobs);
xbuff = zinst*xbuff;
theta = inv(xbuff'*xbuff)*(xbuff'*yobs);
theta = inv(xobs'*xobs)*(xobs'*yobs);
theta;

__con = 0 ;
__altnam = "RC"|"MC"|"SIGMA"|"mtilde" ;
{vnam,m,best,stb,vc,stderr,sigma,cx,rsq,resid,dwstat} = ols(0,yobs,xobs);

end ;