JointMLEMathProgEC.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/jointMC.dat";

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

# logging options
let debug_log := 1;

# Define the MLE problem
problem JointMarkovActvMLE:

# Choose the objective function
likelihood_joint,

# List the variables
EW, 
actvUtil,
sumActvUtil,
jointChoiceUtil,
jointChoiceProb,
Um, b, c, rho,
# EV, actvUtil, choiceUtil, choiceProb, VoT, theta, Uw, xi, gamma,

# List the constraints
# Symmetric_Um,
# Symmetric_b,
# Symmetric_c,
Bellman_Joint,
Bellman_JointH;

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

# Set at a trivial initial value
let {t in 0..H, (j1,j2) in AW1xAW2} EW[t,j1,j2] := 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;

# let {n in PERS, j in ACTV} Um[n,j] := Um0[n,j];
# let {n in PERS, j in ACTV}  b[n,j] := b0[n,j];
# fix {n in PERS, j in ACTV}  c[n,j] := c0[n,j];
# let {j in ALLACTV} rho[j] := rho0[j];

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

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

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

# export to MATLAB
printf "ml = %f;\n", likelihood_joint > DATA/jointMLE.m;

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

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

# display the value function
# display EV;