// ------------------------------------------------------------------------------ // This file solves the dynamic model in // "Estimating Counterfactuals in Dynamic Discrete Choice // Structural Models Using Euler Equations" // using five different solution methods // // Tables 2 and 3 in the paper displayed at end of program for // (user) selected state space and persistence of // transition probabilities // // Victor Aguirregabiria and Arvind Magesan // November 2016 // ------------------------------------------------------------------------------ new; library pgraph; seed =2393457892; //--------------------------------------------------------------------------------// // Contents // //--------------------------------------------------------------------------------// //0. Constants and primitives //1. Transition of (exogenous) state variables //2. The state space //3. Transition of (endogenous) state variable //4. Function for Valuation step in Policy iterations (Newton-Kantorovich) //--------------------------Monte Carlo Simulations Start Here------------------------- //5. Solving the model using PF iterations (STPI) //6. Solving the model using EE-P iterations (EEPI) //7. Solving the model using Value iterations (STVF) //8. Solving the model using EE-V iterations (EEVF) //9. Solving the model using Relative-Value iterations (RVF) //----------------------------Monte Carlo Simulations End Here------------------------- //Set which solution methods are "switched on." If value set to zero the program will not solve using that method. STPI = 1; //Policy Iteration (Newton-Kantorovich) EEPI = 1; //Euler Equations Policy Iteration STVF = 1; //Standard Value Iteration EEVF = 1; //Euler Equations Value Iteration RVF = 1; //Relative Value Iteration //zsize sets the dimensionality of the problem, dimx: |X| = 2*zsize^5 zsize = 3; dimx = 2*zsize^(5); //sigek sets the persistence of the transition probabilities (explained below). Lower sigek means more persistent, sigek must be greater than 0. sigek = 1; /****************/ /* 0. Constants */ /****************/ beta = 0.95 ; @ Discount factor @ //The payoff to being active is given by: //\pi(1, \mathbf{x}_t) = {VP}_{it} - {FC}_{it} - {EC}_{it}*(1-y_{it}) //where: // @ Variable Profit Specification @ VP_{it} = p_t \exp(\lambda \omega_{it} ) // @ Price Specification @ p_t = \alpha_0^{VP} + \alpha_1^{VP} z_{1t} + \alpha_2^{VP} z_{2t} // @ Fixed Cost Parameter @ FC_{it} = \alpha_0^{FC} + \alpha_1^{FC}z_{3t} // @ Entry Cost Parameter @ EC_{it} = \alpha_0^{EC} + \alpha_1^{EC}z_{4t} lamda = 1; alpha0_VP = 0.5; alpha1_VP = 1; alpha2_VP = -1; alpha0_FC = 1.5; alpha1_FC = 1; alpha0_EC = 1; alpha1_EC = 1; sigmae = 1; @ Structural error parameter@ theta_0 = lamda|alpha0_VP|alpha1_VP|alpha2_VP|alpha0_FC|alpha1_FC|alpha0_EC|alpha1_EC|sigmae; @ Payoff Parameter Vector@ //Parameter vector in identified form theta_ID = (alpha0_VP/sigmae)|(alpha1_VP/sigmae)|(alpha2_VP/sigmae)|(alpha0_FC/sigmae)|(alpha1_FC/sigmae)|(alpha0_EC/sigmae)|(alpha1_EC/sigmae); tol = 1e-8; /*************************************/ /* 1. Transition of State Variables */ /*************************************/ //We assume the exogenous variables z1, z2, z3, z4, omega follow independent AR(1) processes // z_kt = \gammak_0 + \gammak_1 z_{kt-1} + e_kt // e_kt ~ N(0,sigek) //The parameters sigek control the persistence of the transition probability matrix. Higher values imply more variability in Z over time, and thus less persistence. //This function takes as inputs parameters of an AR(1) process //and the grid and constructs the matrix of transition probabilities // ----------------------------------------------------------------- // // Format {transition} = tranfun(grid, par) // // Input grid - vector of observations of the dependent variable // par - vector of parameters including error variance. Intercept, slope, error variance. // // Output transition - dim(grid) X dim(grid) matrix of transition probabilities // ----------------------------------------------------------------- // proc (1) = tranfun(grid, par) ; local w, tran, numv, r0, c0, gam0, gam1, sige ; numv = rows(grid); w = (maxc(grid)-minc(grid))/(numv-1); tran = zeros(numv, numv); gam0 = par[1]; gam1 = par[2]; sige = par[3]; r0 = 1; do while r0<=numv; //Loop in todays state c0 = 1; do while c0<=numv; //Loop in tomorrow's state tran[r0,c0] = cdfn((grid[c0]-gam0-gam1*grid[r0]+(w/2))/sqrt(sige) ) - cdfn((grid[c0]-gam0-gam1*grid[r0]-(w/2))/sqrt(sige) ); if c0==1; tran[r0,c0] = cdfn((grid[1]-gam0-gam1*grid[r0]+(w/2))/sqrt(sige) ); endif; if c0==numv; tran[r0,c0] = 1 - cdfn((grid[numv]-gam0-gam1*grid[r0]-(w/2))/sqrt(sige) ); endif; c0 = c0+1; endo; r0 = r0+1; endo; retp(tran) ; endp ; // This function takes as inputs the transition of the five exogenous variables (z_1, z_2, z_3, z_4, omega) // at a given state and returns the joint distribution (product) // ----------------------------------------------------------------- // // Format {tran_out} = tranprod(f1, f2, f3, f4, f5) // // Input f1-f5 - row vector of transition probabilities // Output tran_out - joint probability distribution // ----------------------------------------------------------------- // proc (1) = tranprod(f1, f2, f3, f4, f5) ; local tran_out ; tran_out = (((f1.*.f2).*.f3).*.f4).*.f5; retp(tran_out) ; endp ; //Z1: Market size //Space of Z1 numz1 = zsize; maxz1 = 1; minz1 = 0; w1 = (maxz1-minz1)/(numz1-1); z1val = seqa(minz1,w1, numz1); z1id = seqa(1,1,numz1); sige1 = sigek; gam10 = 0; gam11 = 0.6; par1 = gam10|gam11|sige1; //Call function to calculate Z1 transition matrix {z1tran} = tranfun(z1val, par1); //Z2: Competitors //Space of Z2 numz2 = zsize; maxz2 = 1; minz2 = 0; w2 = (maxz2-minz2)/(numz2-1); z2val = seqa(minz2,w2, numz2); z2id = seqa(1,1,numz2); sige2 = sigek; gam20 = 0; gam21 = 0.6; par2 = gam20|gam21|sige2; //Call function to calculate Z2 transition matrix {z2tran} = tranfun(z2val, par2); //Z3: Fixed Cost Shock //Space of Z3 numz3 = zsize; maxz3 = 1; minz3 = -maxz3; w3 = (maxz3-minz3)/(numz3-1); z3val = seqa(minz3,w3, numz3); z3id = seqa(1,1,numz3); sige3 = sigek; gam30 = 0; gam31 = 0.6; par3 = gam30|gam31|sige3; //Call function to calculate Z3 transition matrix {z3tran} = tranfun(z3val, par3); //Z4: Fixed Cost Shock //Space of Z4 numz4 = zsize; maxz4 = 1; minz4 = -maxz4; w4 = (maxz4-minz4)/(numz4-1); z4val = seqa(minz4,w4, numz4); z4id = seqa(1,1,numz4); sige4 = sigek; gam40 = 0; gam41 = 0.6; par4 = gam40|gam41|sige4; //Call function to calculate Z4 transition matrix {z4tran} = tranfun(z4val, par4); //Productivity Omega //Space of omega numom = zsize; maxom = 1; minom = -1; wom = (maxom-minom)/(numom-1); omval = seqa(minom,wom, numom); omid = seqa(1,1,numom); sigeom = sigek; gamom0 = 0.2; gamom1 = 0.9; parom = gamom0|gamom1|sigeom; //Call function to calculate omega transition matrix {omtran} = tranfun(omval, parom); /*****************************************/ ; /* 2. State Space of Exogenous Variables */ ; /*****************************************/ ; //Create the space and transition of exogenous variables - this is a large matrix zstates = (z1val.*.ones(numz2,1))~(ones(numz1,1).*.z2val); zstates = (zstates.*.ones(numz3,1))~(ones(rows(zstates),1).*.z3val); zstates = (zstates.*.ones(numz4,1))~(ones(rows(zstates),1).*.z4val); numz = rows(zstates); zval = zstates; clear zstates; //Add Omega zval = (zval.*.ones(numom,1))~(ones(numz,1).*.omval); numz = rows(zval); //Construct same matrix for state ID zid = (z1id.*.ones(numz2,1))~(ones(numz1,1).*.z2id); zid = (zid.*.ones(numz3,1))~(ones(rows(zid),1).*.z3id); zid = (zid.*.ones(numz4,1))~(ones(rows(zid),1).*.z4id); zid = (zid.*.ones(numom,1))~(ones(rows(zid),1).*.omid); zid = zid|zid; /***********************************/ ; /* 3. Transition incumbency status */ ; /***********************************/ ; astate = 0|1; atrans0 = (1~0)|(1~0); atrans1 = (0~1)|(0~1); aid = seqa(1,1,rows(astate)); //Add incumbency to get full state space. states = (astate.*.ones(numz,1))~(ones(2,1).*.zval); nstates = rows(states); //States matrix: rows are states, columns are: y~z1~z2~z3~z4~omega //NOTE: WE DO NOT CALCULATE FULL TRANSITION MTX HERE - TO SAVE MEMORY WE CALCULATE AS WE NEED WHILE SOLVING THE MODEL /****************************************************************************/ /* 4. Function for Valuation step in Policy iterations (Newton-Kantorovich) */ /****************************************************************************/ // When This function is called it ouputs the infinite period forward values associated with a policiy // ----------------------------------------------------------------- // // Format { ztilde1, ztilde0, etilde1, etilde0} = STmapping(pay, contran, tran, pinit) // // Input pay - variables in payoffs as they appear in ST mapping, dimension: S x (r). This is "Z1." // contran - An indicator: if set to 1, transition is constructed by // calling another function, if set to 0 we include the transition matrix // tran - If contran = 1 we set tran = 0, if contran = 0 we include the transition matrix (of exogenous variables) here // This choice depends on dimensionality |X| // // Output ztilde1, ztilde0 - choice specific values associated with observable components // etilde1, etilde0 - choice specific values associated with unobservable components // ----------------------------------------------------------------- // proc (4) = STmapping(pay, contran, tran, pinit) ; local numcp, CP_in, it_s, e0, e1, z1, z0, zbar_p, ebar_p, W_Pz0, W_Pe0, W_Pz1, W_Pe1, ztilde0, ztilde1, x0id, exx0tran, etilde0, etilde1, Wit, Wcrit, x0, F1x0, F0x0, Fbar_Px0, F0, F1, Fbar_P, W_Pz, W_Pe; numcp = rows(pay); CP_in = pinit; it_s = 0; do while it_s<1; e0 = 0.5772156649 - ln(1-CP_in); e1 = 0.5772156649 - ln(CP_in); z1 = pay; z0 = zeros(numcp,cols(z1)); zbar_P = CP_in.*z1 + (1-CP_in).*z0; ebar_P = CP_in.*e1 + (1-CP_in).*e0; if contran==1; W_Pz0 = zeros(numcp, cols(z1)); W_Pe0 = zeros(numcp, 1); W_Pz1 = ones(numcp, cols(z1)); W_Pe1 = ones(numcp, 1); ztilde0 = zeros(numcp, cols(z1)); ztilde1 = zeros(numcp, cols(z1)); etilde0 = zeros(numcp, 1); etilde1 = zeros(numcp, 1); Wit=0; Wcrit = maxc(maxc(abs(W_Pz1-W_Pz0)')|abs(W_Pe1-W_Pe0)) ; do while (Wcrit>=1e-8); //Construct the W matrices and then ztilde/etilde matrices state by state to avoid storing transition in memory //WE DON'T NEED TO CALCULATE ztildes and etildes everytime, we do it here to avoid calculating transition again below. We are calculating transition //here anyways. We use the previous iteration's WPZ0 WPE0 to construct ztilde and etilde, at convergence of W's we have the right values. x0 = 1; do while x0<=nstates; x0id = zid[x0,.]; //We are using an object defined globally here. It is only for solving on the full state space {exx0tran} = tranprod(z1tran[x0id[1],.], z2tran[x0id[2],.], z3tran[x0id[3],.], z4tran[x0id[4],.], omtran[x0id[5],.]) ; F1x0 = atrans1[1,.].*.exx0tran; F0x0 = atrans0[1,.].*.exx0tran; Fbar_Px0 = CP_in[x0].*F1x0 + (1-CP_in[x0]).*F0x0; W_Pz1[x0,.] = zbar_P[x0,.] + beta*Fbar_Px0*W_Pz0; W_Pe1[x0,.] = ebar_P[x0,.] + beta*Fbar_Px0*W_Pe0; //Calculate ztilde and etildes only when close to convergence if Wcrit<=1; ztilde0[x0,.] = z0[x0,.] + beta*F0x0*W_Pz0; ztilde1[x0,.] = z1[x0,.] + beta*F1x0*W_Pz0; etilde0[x0,.] = beta*F0x0*W_Pe0; etilde1[x0,.] = beta*F1x0*W_Pe0; endif; x0 = x0+1; endo; Wcrit = maxc(maxc(abs(W_Pz1-W_Pz0)')|abs(W_Pe1-W_Pe0)) ; W_Pz0 = W_Pz1; W_Pe0 = W_Pe1; Wit = Wit+1; endo; else; F0 = atrans0.*.tran; F1 = atrans1.*.tran; Fbar_P = CP_in.*F1 + (1-CP_in).*F0; //Solving W_Pz and W_Pe by method of successive approximation W_Pz0 = zeros(numcp, cols(z1)); W_Pe0 = zeros(numcp, 1); W_Pz1 = ones(numcp, cols(z1)); W_Pe1 = ones(numcp, 1); Wit=0; Wcrit = maxc(maxc(abs(W_Pz1-W_Pz0)')|abs(W_Pe1-W_Pe0)) ; do while (Wcrit>=1e-8); W_Pz1 = zbar_P + beta*Fbar_P*W_Pz0; W_Pe1 = ebar_P + beta*Fbar_P*W_Pe0; Wcrit = maxc(maxc(abs(W_Pz1-W_Pz0)')|abs(W_Pe1-W_Pe0)) ; W_Pz0 = W_Pz1; W_Pe0 = W_Pe1; Wit = Wit+1; endo; W_Pz = W_Pz0; W_Pe = W_Pe0; ztilde0 = z0 + beta*F0*W_Pz; ztilde1 = z1 + beta*F1*W_Pz; etilde0 = beta*F0*W_Pe; etilde1 = beta*F1*W_Pe; endif; it_s = it_s+1; endo; retp(ztilde1, ztilde0, etilde1, etilde0) ; endp ; //------------------------------------------------Monte Carlo Simulations Start Here-------------------------------------------------------------------// nvec=10; initmtx = rndus(nstates,nvec,seed); //Different initial CP vectors statistics = {}; //Collect here statistics on convergence properties for each solution method /**********************************************/ ; /* 5. Compute CPS with standard iterations */ ; /**********************************************/ ; if STPI==1; // For each vector find the fixed point using the standard mapping standard_FP = {}; standard_iterations = {}; standard_time = {}; standard_const = {}; v0 = 1; do while v0<=nvec; CP_0 = initmtx[.,v0]; CP_1 = ones(nstates,1); standard_CP = initmtx[.,v0]; crit = maxc(abs(CP_1-CP_0)) ; crit_s = crit; it_s = 0; //Flow Payoff components : Variable Profit, fixed cost, entry cost price = ones(nstates,1)~states[.,2]~states[.,3]; VP = price.*exp(lamda*states[.,6] ); FC = ones(nstates,1)~states[.,4]; EC = ones(nstates,1)~states[.,5]; z1 = VP~(-FC)~(-EC.*(1-states[.,1])); z0 = zeros(nstates,cols(z1)); timein = {}; do while (crit>=tol); time0in = time; { ztilde1, ztilde0, etilde1, etilde0} = STmapping(z1, 1, 0, CP_0); ztilde0u = (ztilde0)*theta_ID ; ztilde1u = (ztilde1)*theta_ID ; kern = ztilde1u-ztilde0u + etilde1-etilde0 ; CP_1 = exp(kern)./(1+exp(kern)); time1in = time; timein = timein|((time1in-time0in)'*(3600 | 60 | 1| 0.01)); standard_CP = standard_CP~CP_1; crit = maxc(abs(CP_1-CP_0)) ; crit_s = crit_s|crit; CP_0 = CP_1; it_s = it_s+1; endo; time_s =meanc(timein); standard_FP = standard_FP~CP_0; standard_iterations = standard_iterations|it_s; standard_time = standard_time|time_s; constk0 = {}; k0 = 2; do while k0<=cols(standard_CP)-1; constk0 = constk0|(maxc(abs(standard_CP[.,k0+1] - standard_CP[.,k0]))./maxc(abs(standard_CP[.,k0] - standard_CP[.,k0-1]))); k0 = k0+1; endo; standard_const = standard_const|maxc(constk0); v0 = v0+1; endo; statistics = statistics|(meanc(standard_iterations)~meanc(standard_time)~(meanc(standard_time.*standard_iterations))~maxc(standard_const)); endif; /***********************************************/ ; /* 6. Compute CPS using EE-P iterations */ ; /***********************************************/ ; if EEPI==1; EE_FP = {}; EE_iterations = {}; EE_time = {}; EE_const = {}; v0 = 1; do while v0<=nvec; CPcf_0 = initmtx[.,v0]; CPcf_1 = ones(nstates,1); EE_CP = initmtx[.,v0]; crit = maxc(abs(CPcf_1-CPcf_0)) ; crit_e = {}; it_e = 0; price = ones(nstates,1)~states[.,2]~states[.,3]; VP = price.*exp(lamda*states[.,6] ); FC = ones(nstates,1)~states[.,4]; EC = ones(nstates,1)~states[.,5]; timein = {}; do while (crit>=tol); time0in = time; log_diff = ln(1-CPcf_0[nstates/2+1:nstates]) - ln(1-CPcf_0[1:nstates/2]); //For calculating expectations of future variables we go state by state to avoid storing a large matrix //Since we are integrating over exogenous objects only we only need to loop through exogenous states elog_diff = {}; x0 = 1; do while x0<=numz; x0id = zid[x0,.]; fz = (((z1tran[x0id[1],.].*.z2tran[x0id[2],.]).*.z3tran[x0id[3],.]).*.z4tran[x0id[4],.]).*.omtran[x0id[5],.]; elog_diff = elog_diff|(fz*log_diff); x0 = x0+1; endo; elog_diff = elog_diff|elog_diff; cfkern = (VP~(-FC)~(-EC.*(1-states[.,1])))*theta_ID-beta*elog_diff ; CPcf_1= exp(cfkern)./(1+exp(cfkern)); time1in = time; timein = timein|(time1in-time0in)'*(3600 | 60 | 1| 0.01); time_e = meanc(timein); EE_CP = EE_CP~CPcf_1; crit = maxc(abs(CPcf_1-CPcf_0)) ; crit_e = crit_e|crit; CPcf_0 = CPcf_1; it_e = it_e+1; endo; EE_FP = EE_FP~CPcf_0; EE_iterations = EE_iterations|it_e; EE_time = EE_time|time_e; constk0 = {}; k0 = 2; do while k0<=cols(EE_CP)-1; constk0 = constk0|(maxc(abs(EE_CP[.,k0+1] - EE_CP[.,k0]))./maxc(abs(EE_CP[.,k0] - EE_CP[.,k0-1]))); k0 = k0+1; endo; EE_const = EE_const|maxc(constk0); v0 = v0+1; endo; statistics = statistics|(meanc(EE_iterations)~meanc(EE_time)~(meanc(EE_time.*EE_iterations))~maxc(EE_const)); endif; /*********************************************/ ; /* 7. Compute CPS using value iteration */ ; /*********************************************/ ; if STVF==1; initmtxVF = rndus(nstates,nvec,seed); VF_FP = {}; //Keep the fixed points associated with each initial vector VF_iterations = {}; VF_time = {}; VF_const = {}; v0 = 1; do while v0<=nvec; VF_0 = initmtxVF[.,v0]; VF_1 = ones(nstates,1); VFmtx = VF_0; crit = maxc(abs(VF_1-VF_0)) ; crit_VF = crit; it_VF = 0; price = ones(nstates,1)~states[.,2]~states[.,3]; VP = price.*exp(lamda*states[.,6] ); FC = ones(nstates,1)~states[.,4]; EC = ones(nstates,1)~states[.,5]; z1 = VP~(-FC)~(-EC.*(1-states[.,1])); z0 = zeros(nstates,cols(z1)); pay1 = z1*theta_ID; pay0 = zeros(rows(pay1),cols(pay1)); timein = {}; do while (crit_VF>=tol); value0_final = {}; value1_final = {}; numz = nstates/2; x0 = 1; time0in = time; do while x0<=numz; x0id = zid[x0,.]; fz = (((z1tran[x0id[1],.].*.z2tran[x0id[2],.]).*.z3tran[x0id[3],.]).*.z4tran[x0id[4],.]).*.omtran[x0id[5],.]; fzv0 = beta * fz * VF_0[1:numz]; fzv1 = beta * fz * VF_0[numz+1:nstates]; CVF00_x0 = pay0[x0] + fzv0; CVF10_x0 = pay1[x0] + fzv1; CVF01_x0 = pay0[x0+numz] + fzv0; CVF11_x0 = pay1[x0+numz] + fzv1; value0_final = value0_final|ln(exp(CVF00_x0) + exp(CVF10_x0) ); value1_final = value1_final|ln(exp(CVF01_x0) + exp(CVF11_x0) ); x0 = x0+1; endo; time1in = time; timein = timein|((time1in-time0in)'*(3600 | 60 | 1| 0.01)); VF_1 = value0_final|value1_final; VFmtx = VFmtx~VF_1; crit = maxc(abs(VF_1-VF_0)) ; crit_VF = crit_VF|crit; VF_0 = VF_1; it_VF = it_VF+1; endo; time_VF = meanc(timein); VF_iterations = VF_iterations|it_VF; VF_time = VF_time|time_VF; constk0 = {}; k0 = 2; do while k0<=cols(VFmtx)-1; constk0 = constk0|(maxc(abs(VFmtx[.,k0+1] - VFmtx[.,k0]))./maxc(abs(VFmtx[.,k0] - VFmtx[.,k0-1]))); k0 = k0+1; endo; VF_const = VF_const|maxc(constk0); v0 = v0+1; endo; statistics = statistics|(meanc(VF_iterations)~meanc(VF_time)~(meanc(VF_time.*VF_iterations))~maxc(VF_const)); endif; /****************************************************/ ; /* 8. Compute CPS using EE - Value mapping */ ; /****************************************************/ ; if EEVF==1; EEV_FP = {}; //Keep the fixed points associated with each initial vector EEV_iterations = {}; EEV_time = {}; EEV_const = {}; v0 = 1; do while v0<=nvec; EEv_0 = initmtx[.,v0]; EEv_1 = ones(nstates,1); EEVF = initmtx[.,v0]; //Store the steps along the way to convergence to calculate Lipschitz constant crit = maxc(abs(EEv_1-EEv_0)) ; crit_ev = {}; it_ev = 0; timein = {}; price = ones(nstates,1)~states[.,2]~states[.,3]; VP = price.*exp(lamda*states[.,6] ); FC = ones(nstates,1)~states[.,4]; EC = ones(nstates,1)~states[.,5]; pitilde = ((VP~(-FC)~(-EC.*(1-states[.,1]))))*theta_ID; do while (crit>=tol); //We calculate expectations of future variables state by state to avoid storing a large matrix //Since we are integrating over exogenous objects only we only need to loop through exogenous states time0in = time; log_diff = ln(1+exp(EEv_0[numz+1:nstates])) - ln(1+exp(EEv_0[1:numz])); elog_diff = {}; x0 = 1; do while x0<=numz; x0id = zid[x0,.]; fz = (((z1tran[x0id[1],.].*.z2tran[x0id[2],.]).*.z3tran[x0id[3],.]).*.z4tran[x0id[4],.]).*.omtran[x0id[5],.]; elog_diff = elog_diff|(fz*log_diff); x0 = x0+1; endo; elog_diff = elog_diff|elog_diff; EEv_1 = pitilde + beta*elog_diff; time1in = time; timein = timein|(time1in-time0in)'*(3600 | 60 | 1| 0.01); EEVF = EEVF~EEv_1; crit = maxc(abs(EEv_1-EEv_0)) ; crit_ev = crit_ev|crit; EEv_0 = EEv_1; it_ev = it_ev+1; endo; //Solve for CPs CPEE_v = exp(EEv_0)./(1+exp(EEv_0)); EEV_FP = EEV_FP~CPEE_v; EEV_iterations = EEV_iterations|it_ev; time_ev = meanc(timein); EEV_time = EEV_time|time_ev; constk0 = {}; k0 = 2; do while k0<=cols(EEVF)-1; constk0 = constk0|(maxc(abs(EEVF[.,k0+1] - EEVF[.,k0]))./maxc(abs(EEVF[.,k0] - EEVF[.,k0-1]))); k0 = k0+1; endo; EEV_const = EEV_const|maxc(constk0); v0 = v0+1; endo; statistics = statistics|(meanc(EEV_iterations)~meanc(EEV_time)~(meanc(EEV_time.*EEV_iterations))~maxc(EEV_const)); endif; /****************************************************/ ; /* 9. Compute CPS using Relative Value mapping */ ; /****************************************************/ ; if RVF==1; initmtxRVF = rndus(nstates,nvec,seed); RVF_FP = {}; //Keep the fixed points associated with each initial vector RVF_iterations = {}; RVF_time = {}; RVF_const = {}; v0 = 1; do while v0<=nvec; RVF_0 = initmtxRVF[.,v0]; RVF_1 = ones(nstates,1); RVFmtx = RVF_0; crit = maxc(abs(RVF_1-RVF_0)) ; crit_RVF = crit; it_RVF = 0; timein = {}; //Payoff components price = ones(nstates,1)~states[.,2]~states[.,3]; VP = price.*exp(lamda*states[.,6] ); FC = ones(nstates,1)~states[.,4]; EC = ones(nstates,1)~states[.,5]; z1 = VP~(-FC)~(-EC.*(1-states[.,1])); z0 = zeros(nstates,cols(z1)); pay1 = z1*theta_ID; pay0 = zeros(nstates,1); do while (crit_RVF>=tol); //Construct vectors of choice-specific values state by state to avoid storing transition in memory time0in = time; value0_final = {}; value1_final = {}; x0 = 1; do while x0<=numz; x0id = zid[x0,.]; fz = (((z1tran[x0id[1],.].*.z2tran[x0id[2],.]).*.z3tran[x0id[3],.]).*.z4tran[x0id[4],.]).*.omtran[x0id[5],.]; fzv0 = beta * fz * RVF_0[1:numz]; fzv1 = beta * fz * RVF_0[numz+1:nstates]; CVF00_x0 = pay0[x0] + fzv0; CVF10_x0 = pay1[x0] + fzv1; CVF01_x0 = pay0[x0+numz] + fzv0; CVF11_x0 = pay1[x0+numz] + fzv1; value0_final = value0_final|ln(exp(CVF00_x0) + exp(CVF10_x0) ); value1_final = value1_final|ln(exp(CVF01_x0) + exp(CVF11_x0) ); x0 = x0+1; endo; VF_1 = value0_final|value1_final; RVF_1 = VF_1 - VF_1[1].*ones(rows(VF_1),1); time1in = time; timein = timein|(time1in-time0in)'*(3600 | 60 | 1| 0.01); RVFmtx = RVFmtx~RVF_1; crit = maxc(abs(RVF_1-RVF_0)) ; crit_RVF = crit_RVF|crit; RVF_0 = RVF_1; it_RVF = it_RVF+1; endo; time_RVF = meanc(timein); RVF_iterations = RVF_iterations|it_RVF; RVF_time = RVF_time|time_RVF; constk0 = {}; k0 = 2; do while k0<=cols(RVFmtx)-1; constk0 = constk0|(maxc(abs(RVFmtx[.,k0+1] - RVFmtx[.,k0]))./maxc(abs(RVFmtx[.,k0] - RVFmtx[.,k0-1]))); k0 = k0+1; endo; RVF_const = RVF_const|maxc(constk0); v0 = v0+1; endo; statistics = statistics|(meanc(RVF_iterations)~meanc(RVF_time)~(meanc(RVF_time.*RVF_iterations))~maxc(RVF_const)); endif; //Table 2 in paper for state space dimx ""; "Lipschitz Constants"; ""; " EE-v EE-p PF VF RVF"; ""; statistics[4,4]~statistics[2,4]~statistics[1,4]~statistics[3,4]~statistics[5,4]; //Table 3 in paper for state space dimx ""; "Number of Iterations"; ""; " EE-v EE-p PF VF RVF"; ""; statistics[4,1]~statistics[2,1]~statistics[1,1]~statistics[3,1]~statistics[5,1]; ""; "Time Per Iteration"; ""; " EE-v EE-p PF VF RVF"; ""; statistics[4,2]~statistics[2,2]~statistics[1,2]~statistics[3,2]~statistics[5,2]; ""; "Total Time"; ""; " EE-v EE-p PF VF RVF"; ""; statistics[4,3]~statistics[2,3]~statistics[1,3]~statistics[3,3]~statistics[5,3]; ""; "Total Time"; ""; " EE-v EE-p PF VF RVF"; ""; statistics[4,3]~statistics[2,3]~statistics[1,3]~statistics[3,3]~statistics[5,3]; ""; "Time Ratio"; ""; " EE-v/EE-v EE-p/EE-v PF/EE-v VF/EE-v RVF/EE-v"; ""; statistics[4,3]/statistics[4,3]~statistics[2,3]/statistics[4,3]~statistics[1,3]/statistics[4,3]~statistics[3,3]/statistics[4,3]~statistics[5,3]/statistics[4,3]; end;