#####  HAROLD ZURCHER BUS REPAIR EXAMPLE  #####
#
# An NLP method to compute maximum likelihood estimates of the 
# Harold Zurcher bus problem in Rust, Econometrica, 1987. 
# K. Judd and C.-L. Su, 
# February 21, 2006. Current version: October 26, 2006
#
#
#####  SET UP THE PROBLEM  #####
#
#  Define the state space used in the dynamic programming part
param    N; 	# number of states used in dynamic programming approximation
set X := 1..N; 	# X is the index set of states
#  x[i] denotes state i; the set of states is a uniform grid on the interval [xmin, xmax]
param xmin := 0;
param xmax := 100;
param x {i in X} := xmin + (xmax-xmin)/(N-1)*(i-1); 
#
#  Define and process the data
param nT;		# number of periods in data
set T := 1..nT;	# T is the vector of time indices
param Xt {T}; 	# Xt[t] is the true mileage at time t
param dt {T};	# decision at time t
param  Xtmax := max {t in T} Xt[t];    # scale parameter

# scale the observed mileage to the range of [0, 1] for the purpose of scaling the data moments 
param Xts {t in T} := Xt[t]/Xtmax;     
#  The dynamic programming model in the estimation lives on a discrete state
#  Binning process: assign true mileage Xt[t] to the closest state in X
param xt {t in T} := ceil(Xt[t]/(xmax-xmin)*(N-1)+0.5); 

# scale the state x[i] by the same data moment scaling parameter  
param xs {i in X} := x[i]/Xtmax;      
#
#  Define "known" structural parameters
#  We fix beta since data cannot identify it
param beta;      	# discount factor
#
#####  DECLARE STRUCTURAL PARAMETERS TO BE ESTIMATED #####
#####  Parameters for cost function  #####
#  c(x, thetaCost) = thetaCost[1]*x + thetaCost[2]*x^2
var thetaCost {1..2} >= 0; 
#
#####  Parameters and definition of transition process  #####
#  thetaProbs defines Markov chain
var thetaProbs {1..3} >= 0.00001; 
#
#  Define the Markov chain representing the changes in mileage on the x[i] grid.
#
#  The state increases by some amount in [0,JumpMax] where 
#  JumpMax is the maximum increase in mileage in one period
#  and equals a fraction JumpRatio of the range of mileage in the state space
param JumpRatio;             
param JumpMax := (xmax-xmin) * JumpRatio;
#
#  We assume that the jumps are independent of the current state
#  and its transition process has three pieces, each a uniform distribution    
#  Define 1st break point for stepwise uniform distribution in mileage increase
param M1 := ceil(1/4*JumpMax/(xmax-xmin)*(N-1)+0.5);   
#  Define 2nd break point for stepwise uniform distribution in mileage increase
param M2 := ceil(3/4*JumpMax/(xmax-xmin)*(N-1)+0.5);   
#  Define end point for stepwise uniform distribution in mileage increase
param M := ceil(JumpMax/(xmax-xmin)*(N-1)+0.5);        
#
#  Y is the vector of elements in transition rule
set Y := 1..M;
var TransProb {i in Y} = 
	if 	i <= M1 		then 	thetaProbs[1]/M1 
	else if 	i > M1 and i <= M2 	then 	thetaProbs[2]/(M2-M1) 
	else 	thetaProbs[3]/(M-M2);
#
#####  Scrap value parameter  #####
var RC >= 0;
##### END OF STRUCTURAL VARIABLES  #####
#
#####  DECLARE EQUILIBRIUM CONSTRAINT VARIABLES #####
# The NLP approach requires us to solve equilibrium constraint variables
var EV {X};        		# Value Function of each state
#####  END OF EQUILIBRIUM CONSTRAINT VARIABLES  #####
#
#####  DECLARE AUXILIARY VARIABLES  ######
#  Define auxiliary variables to economize on expressions	
#  Create Cost variable to represent the cost function; 
#  Cost[i] is the cost of regular maintenance at x[i].
var Cost {i in X} = sum {j in 1..2} thetaCost[j]*x[i]^(j);	    
#  Let CbEV[i] represent - Cost[i] + beta*EV[i]; 
#  this is the expected payoff at x[i] if regular maintenance is chosen
var CbEV {i in X} = - Cost[i] + beta*EV[i];                     
#  Let PayoffDiff[i] represent -CbEV[i] - RC + CbEV[1]; 
#  this is the difference in expected payoff at x[i] between engine replacement and regular maintenance
var PayoffDiff {i in X} = -CbEV[i] - RC + CbEV[1];               
#  Let ProbRegMaint[i] represent 1/(1+exp(PayoffDiff[i])); 
#  this is the probability of performing regular maintenance at state x[i];
var ProbRegMaint {i in X} = 1/(1+exp(PayoffDiff[i]));     

var BellmanViola {i in 1..(N-M+1)} = sum {j in 0..(M-1)} log(exp(CbEV[i+j])+ exp(-RC + CbEV[1]))* TransProb[j+1] - EV[i];

set D := 1..2;
set S := X cross D;

param d {i in D} := if i = 1 then 0 else 1;

# Define Markov Transition Matrix
var PI {(i1, i2) in S, (j1, j2) in S} = 
    if d[i2] = 0 then 
         if j1 < i1  or j1 - i1 >= M   then    0 
         else if i1 <= N-(M-1) then TransProb[j1-i1+1]*(ProbRegMaint[j1]*(1-d[j2]) +(1-ProbRegMaint[j1])*d[j2]) 
                 else if i1 = N then (ProbRegMaint[N]*(1-d[j2]) +(1-ProbRegMaint[N])*d[j2]) 
                     else if j1 < N then TransProb[j1-i1+1]*(ProbRegMaint[j1]*(1-d[j2]) +(1-ProbRegMaint[j1])*d[j2]) 
                              else (1-sum{k in 0..(N-i1-1)}TransProb[k+1]) *(ProbRegMaint[N]*(1-d[j2]) +(1-ProbRegMaint[N])*d[j2]) 
    else 
         if j1 > M then 0 
         else                    TransProb[j1 -1+1]*(ProbRegMaint[j1]*(1-d[j2]) +(1-ProbRegMaint[j1])*d[j2]) ;

# Define Markov Transition Matrix
# var PI {(xi, di) in S, (xj, dj) in S} = 
#    if di = 0 then 
#         if xj < xi  or xj - xi >= M   then    0 
#         else if xi <= N-(M-1) then TransProb[xj-xi+1]*(ProbRegMaint[xj]*(1-dj) +(1-ProbRegMaint[xj])*dj) 
#                 else if xi = N then (ProbRegMaint[N]*(1-dj) +(1-ProbRegMaint[N])*dj) 
#                     else if xj < N then TransProb[xj-xi+1]*(ProbRegMaint[xj]*(1-dj) +(1-ProbRegMaint[xj])*dj) 
#                              else (1-sum{k in 0..(N-xi-1)}TransProb[k+1]) *(ProbRegMaint[N]*(1-dj) +(1-ProbRegMaint[N])*dj) 
#    else 
#         if xj > M then 0 
#         else                    TransProb[xj -1+1]*(ProbRegMaint[xj]*(1-dj) +(1-ProbRegMaint[xj])*dj) ;

# Define Ergodic distribution vector
var ergoDist {(i1, i2) in S};

# Define moments implied by ergodic distribution 
# mean
var mx = sum {(i, j) in S} ergoDist[i, j]*xs[i];
var md = sum {(i, j) in S} ergoDist[i, j]*d[j];

# variance
var mxx = sum {(i,j) in S} ( ergoDist[i,j]*(xs[i]-mx)^2 );
var mxd = sum {(i,j) in S} ( ergoDist[i,j]*(xs[i]-mx)*(d[j]-md) );
var mdd = sum {(i,j) in S} ( ergoDist[i,j]*(d[j]-md)^2 );

# skewness
var mxxx = sum {(i,j) in S} ( ergoDist[i,j]*(xs[i]-mx)^3 );
var mxxd = sum {(i,j) in S} ( ergoDist[i,j]*(xs[i]-mx)^2*(d[j]-md) );
var mxdd = sum {(i,j) in S} ( ergoDist[i,j]*(xs[i]-mx)*(d[j]-md)^2 );
var mddd = sum {(i,j) in S} ( ergoDist[i,j]*(d[j]-md)^3 );

# serial covariance
# 1 period
# var SeriCOVxx1 = sum {(i1,i2) in S} ( ergoDist[i1,i2]*(xs[i1]-mx) * (sum {(j1, j2) in S} PI[i1, i2, j1, j2] * (xs[j1] - mx) )  );
# var SeriCOVxd1 = sum {(i1,i2) in S} ( ergoDist[i1,i2]*(xs[i1]-mx) * (sum {(j1, j2) in S} PI[i1, i2, j1, j2] * (d[j2] - md) )  );
# var SeriCOVdd1 = sum {(i1,i2) in S} ( ergoDist[i1,i2]*(d[i2]-md) * (sum {(j1, j2) in S} PI[i1, i2, j1, j2] * (d[j2] - md) )  );

# 2 period
# var SeriCOVxx2 = sum {(i1,i2) in S} ( ergoDist[i1,i2]*(xs[i1]-mx) * (sum {(j1, j2) in S, (k1, k2) in S} PI[i1, i2, j1, j2] * PI[j1, j2, k1, k2]*(xs[k1] - mx) )  );
# var SeriCOVxd2 = sum {(i1,i2) in S} ( ergoDist[i1,i2]*(xs[i1]-mx) * (sum {(j1, j2) in S, (k1, k2) in S} PI[i1, i2, j1, j2] * PI[j1, j2, k1, k2]*(d[k2] - md) )  );
# var SeriCOVdd2 = sum {(i1,i2) in S} ( ergoDist[i1,i2]*(d[i2]-md) * (sum {(j1, j2) in S, (k1, k2) in S} PI[i1, i2, j1, j2] * PI[j1, j2, k1, k2]*(d[k2] - md) )  );

# Define data moments
# mean
param MX = (sum {t in T} Xts[t] )/nT;
param MD = (sum {t in T} dt[t] )/nT;

# variance
param MXX =  (sum {t in T}  (Xts[t] - MX)^2)  / nT;
param MXD =  (sum {t in T}  (Xts[t] - MX)*(dt[t] - MD))  / nT;
param MDD =  (sum {t in T}  (dt[t] - MD)^2)  / nT;

# skewness
param MXXX =  (sum {t in T}  (Xts[t] - MX)^3)  / nT;
param MXXD =  (sum {t in T}  (Xts[t] - MX)^2*(dt[t] - MD))  / nT;
param MXDD =  (sum {t in T}  (Xts[t] - MX)*(dt[t] - MD)^2 )  / nT;
param MDDD =  (sum {t in T}  (dt[t] - MD)^3)  / nT;


#####  END OF AUXILIARY VARIABLES  #####
#
#####  OBJECTIVE AND CONSTRAINT DEFINITIONS  #####
#
## Define objective function
##
## # Likelihood function 
# maximize Likelihood: 
## #  The likelihood function contains two pieces
## #  First is the likelihood that the engine is replaced given time t state in the data.
#	sum {t in 2..nT} log(dt[t]*(1-ProbRegMaint[xt[t]]) + (1-dt[t])*ProbRegMaint[xt[t]])
## #  Second is the likelihood that the observed transition between t-1 and t would have occurred.
#	+ sum {t in 2..nT} log(dt[t-1]*(TransProb[xt[t]-1+1]) + (1-dt[t-1])*(TransProb[xt[t]-xt[t-1]+1])); 

minimize GMM: 
#  The moment function 
(MX-mx)^2 + (MD-md)^2 + (MXX-mxx)^2 + (MXD-mxd)^2 + (MDD-mdd)^2 +  (MXXX-mxxx)^2 + (MXXD-mxxd)^2 + (MXDD-mxdd)^2+  (MDDD-mddd)^2 ; 


#
#  List the constraints
#
subject to 
#
#  Bellman equation for states below N-M
	Bellman_1toNminusM {i in X: i <= N-(M-1)}: 
		EV[i] = sum {j in 0..(M-1)} 
			log(exp(CbEV[i+j])+ exp(-RC + CbEV[1]))* TransProb[j+1];
#
#  Bellman equation for states above N-M, (we adjust transition probabilities to keep state in [xmin, xmax])
	Bellman_LastM {i in X: i > N-(M-1) and i <= N-1}: 
		EV[i] = (sum {j in 0..(N-i-1)} 
			log(exp(CbEV[i+j])+ exp(-RC + CbEV[1]))* TransProb[j+1]) 
			+ (1- sum {k in 0..(N-i-1)} TransProb[k+1]) * log(exp(CbEV[N])+ exp(-RC + CbEV[1]));
#
#  Bellman equation for state N
	Bellman_N:  EV[N] = log(exp(CbEV[N])+ exp(-RC + CbEV[1]));
#
#  The probability parameters in transition process must add to one
    Probability: sum {i in 1..3} thetaProbs[i] = 1;

#
#  Put bound on EV; this should not bind, but is a cautionary step to help keep algorithm within bounds
    EVBound {i in X}: EV[i] <= 50;
#

# Constraints for computing Ergodic Distribution

   ergoDistProbability: sum {(i1, i2) in S} ergoDist[i1, i2] = 1;
   ergoEquation {(i1, i2) in S}:   sum{(j1, j2) in S} PI[j1, j2, i1, i2] * ergoDist[j1, j2] = ergoDist[i1, i2];

####   DEFINE THE PROBLEM   #####
#
# Name the problem
problem MPECZurcher:
#
# Choose the objective function
# Likelihood,
GMM,

#
# List the variables
EV, RC, thetaCost, thetaProbs, TransProb, Cost, CbEV, PayoffDiff, ProbRegMaint, BellmanViola, PI, ergoDist,
mx, md, mxx, mxd, mdd, mxxx, mxxd, mxdd, mddd, 

# SeriCOVxx1,  SeriCOVxd1,  SeriCOVdd1,  SeriCOVxx2,  SeriCOVxd2,  SeriCOVdd2, 

#
# List the constraints
Bellman_1toNminusM,
Bellman_LastM, 
Bellman_N, 
Probability,  
EVBound, 
ergoDistProbability,
ergoEquation;
#
#################################