Model code

This part documents the different functions for the calculation of the model objects determining the decision of Harold Zurcher. Following Rust (1987), the code does not estimate the discount factor and it needs to be externally set.

Observed costs

The observed costs are saved in \(num\_states \times 2\) dimensional numpy array. The first column contains the maintenance and the second the replacement costs for each state. The function to calculate the observed costs is:

calc_obs_costs(num_states, maint_func, ...)

Calculating the observed costs of maintenance and replacement for each state.

The inputs are besides the size of the state space, the type of the maintenance cost function as well as the cost parameters and the scale. The inputs will be explained in the following:

Maintenance cost function

So far the code allows for five functional forms. The following table reports the different functional forms for an arbitrary state \(x\). Afterwards I list the APIs of each function and their derivatives. \(states\) is the size of the state space.

Name	Functional form
linear	\(c(x,\theta_1) = \theta_{11} x\)
square root	\(c(x,\theta_1) = \theta_{11} \sqrt{x}\)
cubic	\(c(x,\theta_1) = \theta_{11}x+\theta_{12} x2+\theta_{13} x3\)
hyperbolic	\(c(x,\theta_1) = (\theta_{11} / ((states + 1) - x))\)
quadratic	\(c(x,\theta_1) = (\theta_{11} x +\theta_{12} x**2)\)

Linear cost function x .. currentmodule:: ruspy.model_code.cost_functions

`lin_cost`(num_states, params, scale)	Calculating for each state the observed costs of maintenance in the case of a linear cost function.
`lin_cost_dev`(num_states, scale)	Calculating for each state the derivative of the linear maintenance cost function.

Square root function

`sqrt_costs`(num_states, params, scale)	Calculating for each state the observed costs of maintenance in the case of a square root cost function.
`sqrt_costs_dev`(num_states, scale)	Calculating for each state the derivative of the square root maintenance cost function.

Cubic cost function

`cubic_costs`(num_states, params, scale)	Calculating for each state the observed costs of maintenance in the case of a cubic cost function.
`cubic_costs_dev`(num_states, scale)	Calculating for each state the derivative of the cubic maintenance cost function.

Quadratic cost function

`quadratic_costs`(num_states, params, scale)	Calculating for each state the observed costs of maintenance in the case of a quadratic cost function.
`quadratic_costs_dev`(num_states, scale)	Calculating for each state the derivative of the quadratic maintenance cost function.

hyperbolic cost function

`hyperbolic_costs`(num_states, params, scale)	Calculating for each state the observed costs of maintenance in the case of a hyperbolic cost function.
`hyperbolic_costs_dev`(num_states, scale)	Calculating for each state the derivative of the hyperbolic maintenance cost function.

Cost parameters

The second input are the cost parameters, which are stored as a one dimensional numpy.array. At the first position always the replacement cost \(RC\) is stored. The next positions are subsequently filled with \(\theta_{11}, \theta_{12}, ...\). The exact number depends on the functional form.

Scale

The maintenance costs are, due to feasibility of the fixed point algorithm scaled. The scaling varies across functional forms. The following table contains an overview of the minimal scale needed for each form:

Cost function	Scale
linear	\(10^{-3}\)
square root	\(10^{-2}\)
cubic	\(10^{-8}\)
hyperbolic	\(10^{-1}\)
quadratic	\(10^{-5}\)

Fixed Point Algorithm

This part documents the core contribution to research of the Rust (1987) paper, the Fixed Point Algorithm (FXP). It allows to consequently calculate the log-likelihood value for each cost parameter and thus, to estimate the model and hence builds the corner stone of the Nested Fixed Point Algorithm (NFXP). The computation of the fixed point is managed by:

calc_fixp(trans_mat, obs_costs, disc_fac[, ...])

Calculating the expected value of maintenance fixed point with the polyalgorithm proposed by Rust (1987) and Rust (2000).

The algorithm is set up as a polyalgorithm combining contraction and Newton-Kantorovich (Kantorovich, 1948) iterations. It starts by executing contraction iterations, until it reaches some convergence threshold and then switches to Newton-Kantorovich iterations. The exact mathematical deviation of the separate steps are very nicely illustrated in Rust (2000). The function of these two steps are the following in ruspy:

`contraction_iteration`(ev, trans_mat, ...)	Calculating one iteration of the contraction mapping.
`kantorovich_step`(ev, trans_mat, obs_costs, ...)	Calculating one Newton-Kantorovich step for approximating the fix-point.

Algorithmic details

In the following the variable keys are presented, which allow to specify the algorithmic behavior. The parameters can be grouped into two categories. Switching parameters, which allow to specify, when the algorithm switches from contraction to Newton-Kantorovich iterations and general parameters, which let the algorithm stop. So far, there is no switching back implemented.

max_cont_steps : (int) The maximum number of contraction iterations before switching to Newton-Kantorovich iterations. Default is 20.
switch_tol : (float) If this threshold is reached by contraction iterations, then the algorithm switches to Newton-Kantorovich iterations. Default is \(10^{-3}\).
max_newt_kant_steps : (int) The maximum number of Newton-Kantorovich iterations before the algorithm stops. Default is 20.
threshold : (float) If this threshold is reached by Newton-Kantorovich iterations, then the algorithm stops. Default is \(10^{-12}\).

Expected value of maintenance

In ruspy the expected value of maintenance is stored in a state space sized numpy array. Thus, the exected value of replacement can be found in the zero entry. It is generally denoted by ev, except in the simulation part of the package where it is denoted by ev_known. This illustrates that the expected value is created by the agent on his beliefs of the process.

Common model objects

Here are some common objects with a short description documented.

Transition matrix

The transition matrix for the Markov process are stored in a \(num\_states \times num\_states\) dimensional numpy array. As transition in the case of the replacement corresponds to a transition from state 0, it exists only matrix.

Choice probabilities

The choice probabilities are stored in a \(num\_states \times 2\) dimensional numpy array. So far only choice probabilities, resulting from an unobserved shock with i.i.d. gumbel distribution are implemented. The multinomial logit formula is herefore implemented in:

choice_prob_gumbel(ev, obs_costs, disc_fac)

This function calculates the choice probabilities to maintain or replace the bus engine for each state.

Discount factor

The discount factor, as described in the economic model section, is stored as a float in ruspy. It needs to be set externally for the simulation as well as for the estimation process. The key in the dictionary herefore is always disc_fac.

Demand Function Calculation

The demand function can be derived using the following function get_demand.

get_demand(init_dict, demand_dict, demand_params)

Calculates the implied demand for a range of replacement costs for a certain number of buses over a certain time period.

Based on the estimated structural parameters \(\theta\) obtained from the function estimate which is based on the model specification in Estimation initialization dictionary, one can now derive the implied demand function. In order to do so one has to provide the following inputs:

Initialization Dictionairy

This is the Estimation initialization dictionary needed for the function get_criterion_function. The get_demand function draws the model specifications needed to calculate demand from this.

Demand Dictionary

This dictionairy provides all the necessary information about how the demand function is supposed to look like and how precisely it is supposed to be calculated. It has to hold the following keys:

RC_lower_bound : (float) The lowest replacement cost for which the demand is supposed to be calculated.
RC_upper_bound : (float) The highest replacement cost for which the demand is supposed to be calculated.
demand_evaluations : (int) The grid size of the replacement cost between RC_lower_bound and RC_upper_bound for which the demand level shall be calculated.
tolerance : (float) The stopping tolerance for the fixed point calculation needed to obtain each demand level.
num_periods : (int) Number of months \(T\) for which the expected demand is derived. Consequently, set it to 12 if you want to get the annual expected demand.
num_buses : (int) Number of buses \(M\) for which the demand is calculated.

Demand Parameters

This numpy array contains the structural parameters of the model for which you want to derive the implied demand. This parametrization can come from an estimation procedure with the estimate function but is not limited to that. The first elements are the transition probabilities \(\theta_{30}, \theta_{31}, ...\) and the elements after are \(RC, \theta_{11}, ...\).

Based on those inputs, the get_demand function gives out the following DataFrame.

Demand Results

This pandas DataFrame has the grid of replacement costs defined by RC_lower_bound, RC_upper_bound and demand_evaluations as an index and gives out for each of them in the column demand the expected demand level over the specified time horizon and for the amount of buses in the fleet. It also contains a column success which indicates whether the fixed point algorithm converged successfully.

The use of the get_demand function is shown in the following replication notebook.