// --------------------------------------------------------------------- // 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.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 ;