# Title: Constrained Optimization Approaches for Estimation of Structural Models
#
# Authors: Che-Lin Su and Kenneth L. Judd
# 
#  
#####  HAROLD ZURCHER BUS REPAIR EXAMPLE  #####
#
# A constrained optimization formulation 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
#  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); 
#
#  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];


#####  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])); 

#
#  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;
#

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

#
# List the variables
EV, RC, thetaCost, thetaProbs, TransProb, Cost, CbEV, PayoffDiff, ProbRegMaint, BellmanViola, 

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

#
#################################