I. Importation of the modules


In [1]:
%matplotlib inline
import pandas as pd
import numpy as np
import scipy.stats as stats
import itertools
from math import exp
from IPython.core.pylabtools import figsize
from matplotlib import pyplot as plt
import seaborn as sns
figsize(12.5, 8)
pd.set_option("display.precision", 3)

II. Notations and definitions

1°) Number of states and state transitions

In this paper, the agent observes each bus at a given time $t$, and makes a decision regarding the replacement of the engine or its maintenance as a function of the observed state of the engine (i.e. its mileage).

As in the original paper, we discretize the number of states into intervals of length 5000, and define a transition density which governs the evolution of mileage from one time period to another:

$$ p(x_{t+1} | x_{t} , i_{t} , \theta_{3}) = \left\{ \begin{array}{l l} g(x_{t+1} - x_{t}, \theta_{3}) & \quad \text{if } i_t = 0\\ g(x_{t+1} - 0, \theta_{3}) & \quad \text{if } i_t = 1 \end{array} \right.$$

Here, the function g is defined by a multinomial distribution on the set {0, 1, 2}. In other words, the probability that a bus' mileage will increase by a value x between now and the next maintenance decision is:

  • $\theta_{31}$ if $x \in [0, 5000[$
  • $\theta_{32}$ if $x \in [5000, 10,000[$
  • $\theta_{33} = 1 - (\theta_{31} + \theta_{32})$ if $x \in [10,000, +\infty[$

2°) Decision criteria

The utility function of the decision maker for a single time period is:

$$ u(x_{t}, i , \theta_{1}) = \left\{ \begin{array}{l l} \qquad \enspace -c(x_t, \theta_1) + \epsilon_t(0)& \quad \text{if } i = 0\\ -RC -c(0, \theta_1) + \epsilon_t(1) & \quad \text{if } i = 1 \end{array} \right. \quad \text{(Errors are I.I.D. standard Gumbel)}$$

The agent will chose the decision which maximizes its utility. If the agent is forward looking (i.e. if $\beta \neq 0$), he does not only maximize his utility in the present time period, but also his future utility discounted by $\beta$.

Since $\beta$ cannot be estimated without knowing the utility function, we will later set $\beta$ to its true value.

3°) Cost function

A functional form has to be imposed on the cost function $c(x_t)$.

Here, we assume that the true cost function is a quadratic function of the current state of the bus:

$$ c(x_{t}) = \theta_{11}x + \theta_{12}x^2$$

III. Data-generating process

Bus mileage

The generation of a fake dataset can be done in two different ways:

  • Generating monthly mileage increase from a continuous distribution, discretize the mileage, and recover the state transition parameters ($\theta_3$) using a maximum-likelihood estimation (as done in the paper).
  • Set the vector $\theta_3$ to an arbitrary value, and generate discrete states transitions with corresponding probabilities.

The advantage of the first method is that this is closer to what is actually done in the paper: indeed, the focal goal of the paper is to recover RC and the cost scaling vector of parameters $\theta_1$. If the decision of the agent is based on the discretized state (and not on the actual mileage), then this approach is not noisier than the second one.

For this reason, we will adopt the first approach, and assumes that:

$$ (M_{t+1} - M_{t}) \sim \Gamma(5000, 5000) $$

In [137]:
a = stats.truncnorm(a=0, b=100000, loc=3000, scale = 6000)

lower, upper = 0, 100000
mu, sigma = 5000, 6000
a = stats.truncnorm((lower - mu) / sigma, (upper - mu) / sigma, loc=mu, scale=sigma)

In [139]:
x = np.linspace(0,15000,1000)
plt.plot(x, a.pdf(x))


Out[139]:
[<matplotlib.lines.Line2D at 0x154f2390>]

In [ ]: