MLEMathProgEC.run

# Call solver and give it options

# Load model and data
model "MarkovActv.mod";
data "MarkovActv.dat";

# include the simulated states and choices
data "DATA/MC.dat";

# Is Cauchy distribution used ? 
fix IS_CAUCHY := 1;

# logging options
let debug_log := 1;

# Define the MLE problem
problem MarkovActvMLE:

# Choose the objective function
likelihood,

# List the variables
{t in 0..H, j in AUW[n1]} EV[n1,t,j], 
{t in 0..H, j in AUW[n1]} actvUtil[n1,t,j],
{(t,j) in X[n1], (k,h) in DA[n1,t,j]} sumActvUtil[n1,t,j,k,h],
{(t,j) in X[n1], (k,h) in D[n1,t,j]} choiceUtil[n1,t,j,k,h],
{(t,j) in X[n1], (k,h) in D[n1,t,j]} choiceProb[n1,t,j,k,h],
{j in ACTV} Um[n1,j],
{j in ACTV} b[n1,j],
{j in ACTV} c[n1,j],
# EV, actvUtil, choiceUtil, choiceProb, VoT, theta, Uw, xi, gamma,

# List the constraints
{(t,j) in X[n1]} Bellman_Eqn[n1,t,j],
Bellman_EqnH[n1];

# include the code that define the state and choice set
include MDPStateAction.run;

# Set at a trivial initial value
let {t in 0..H, j in ACTV} EV[n1,t,j] := initEV;

# Fix theta at the true value
# fix theta := trueTheta;

# Initial guesses set at trivial values; probably not good initial guesses
# fix VoT := trueVoT;

# Specify KNITRO solver options:
option solver "C:\Ziena\KNITRO900\knitroampl\knitroampl.exe";
option knitro_options "alg=2 hessopt=1 outlev=3 maxit=600 xtol=0.0000000001 wantsol=1";

# Output commands
option display_round 6, display_width 120;

# fix {j in ACTV} Um[n1,j]:= Um0[n1,j];
# fix {j in ACTV}  b[n1,j]:= b0[n1,j];
# fix {j in ACTV}  c[n1,j]:= c0[n1,j];

# Solve command
solve MarkovActvMLE;
if match (solve_message, "Locally optimal solution") > 0 then {
	printf "solve_message = 'Locally optimal solution';\n" > DATA/MLE.m;
	printf "solve_return = %5.0f;\n", 0 > DATA/MLE.m;
}
else if match (solve_message, "Iteration limit reached") > 0 then {
	printf "solve_message = 'Iteration limit reached';\n" > DATA/MLE.m;
	printf "solve_return = %5.0f;\n", 400 > DATA/MLE.m;
}
else if match (solve_message, "Evaluation error") > 0 then {
	printf "solve_message = 'Evaluation error';\n" > DATA/MLE.m;
	printf "solve_return = %5.0f;\n", 502 > DATA/MLE.m;
}
else {
	printf "solve_message = 'Unknown error';\n", solve_message > DATA/MLE.m;
	printf "solve_return = %5.0f;\n", 1000 > DATA/MLE.m;
}

# display solver running time
# display _solve_time;
# display estimates of structural parameters
display beta, theta, VoT;
display beta, theta, VoT > DATA/MLE.sol;
display Um, b, c;
display Um, b, c > DATA/MLE.sol;
# display the true values
display Um0, b0, c0;
display Um0, b0, c0 > DATA/MLE.sol;
# display Uw, xi, gamma, lambda;

# write the estimates
printf "param Um_[%d,*] := \n", n1 > DATA/MLE.dat;
for {j in ACTV} {
	printf "%d\t%f\n", j, Um[n1,j] > DATA/MLE.dat;
}
printf ";\n" > DATA/MLE.dat;
printf "param  b_[%d,*] := \n", n1 > DATA/MLE.dat;
for {j in ACTV} {
	printf "%d\t%f\n", j, b[n1,j] > DATA/MLE.dat;
}
printf ";\n" > DATA/MLE.dat;
printf "param  c_[%d,*] := \n", n1 > DATA/MLE.dat;
for {j in ACTV} {
	printf "%d\t%f\n", j, c[n1,j] > DATA/MLE.dat;
}
printf ";\n" > DATA/MLE.dat;

# export to MATLAB
printf "ml = %f;\n", likelihood > DATA/MLE.m;
printf "n1 = %d;\n", n1 > DATA/MLE.m;

# export fixed parameters to MLE.m
printf "beta = %f;\n", beta > DATA/MLE.m;
printf "theta = %f;\n", theta > DATA/MLE.m;

# export the estimated parameters to MLE.m
printf "Um_ = [ ", n1 > DATA/MLE.m;
for {j in ACTV} {
	printf "%f ", Um[n1,j] > DATA/MLE.m;
}
printf "];\n" > DATA/MLE.m;
printf "b_ = [ ", n1 > DATA/MLE.m;
for {j in ACTV} {
	printf "%f ", b[n1,j] > DATA/MLE.m;
}
printf "];\n" > DATA/MLE.m;
printf "c_ = [ ", n1 > DATA/MLE.m;
for {j in ACTV} {
	printf "%f ", c[n1,j] > DATA/MLE.m;
}
printf "];\n" > DATA/MLE.m;

# display the value function
# display EV;