The bus engine replacement model of Rust (1987) is the simplest dynamic discrete choice model taken to real data, and the template for every model we estimate later. This class builds it from the Bellman equation of Class 7 to a solver that converges in a handful of Newton steps.
Running the code for this lecture
Every code example below is also a runnable notebook in the course code repository,
in the folder session09-sep22/.
You should have cloned that repository already โ if not, the instructions are in the algorithms and complexity lecture.
Update your copy before the class. Editing a file in place makes git pull refuse to
update it, so discard whatever you changed while experimenting:
cd sb-dse-code
git reset --hard HEAD
git pullNothing in this repository is submitted, so experiment freely โ but copy anything you want to keep out of it first, or commit that to a branch of your own.
Setting up the Python environment is covered in Python.
Homework zeta_zurcher: the bus engine model under the hood
Graded homework: task zeta_zurcher in the class repository. Collect it and work on a
copy in your own repository:
git pull upstream main # collect the task
cp -r tasks/zeta_zurcher solutions/zeta_zurcher # work on the copyTake the model of this class apart: verify its Frรฉchet derivative against finite differences, find the starting points from which NewtonโKantorovich breaks down, and compute the ergodic distribution of mileage under the optimal policy โ the distribution of the data the model generates, which we simulate from on Thursday.
The notebook play_zurcher.ipynb in the task folder carries the model class of this
class, zurcher.py, and the details of the task, step by step. The code of this class,
session09-sep22/zurcher.ipynb in the code repository, is worth having beside you
while you work.
Submit it as a pull request, following the git workflow. The solutions are discussed at the start of the next class.
Rust (1987) and Harold Zurcherยถ
๐ Rust (1987), Econometrica โOptimal Replacement of GMC Bus Engines: An Empirical Model of Harold Zurcherโ is the foundational paper of dynamic structural econometrics. It develops the framework that is used across many fields today, and it does so with five ingredients:
A simple dynamic model: a binary discrete state regenerative optimal stopping problem
A smart solution method: a polyalgorithm with fast NewtonโKantorovich iterations
A coherent econometric specification: unobserved states ensuring a non-degenerate likelihood
A new maximum likelihood estimator: the nested fixed point (NFXP) estimator
A very real application with actual data collected by a real person named Harold Zurcher
The citation count of the paper is still growing, forty years on.
Who cares about Harold Zurcher?
The same model, with different names for the state and the choice, sits underneath
Occupational choice: Keane & Wolpin (1997), JPE
Retirement: Rust & Phelan (1997), Econometrica
Brand choice and advertising: Erdem & Keane (1996), Marketing Science
Choice of college major: Arcidiacono (2004), Journal of Econometrics
Individual migration decisions: Kennan & Walker (2011), Econometrica
High school attendance and work decisions: Eckstein & Wolpin (1999), Econometrica
Sales and dynamics of consumer inventory behavior: Hendel & Nevo (2006), Econometrica
Advertising, learning, and consumer choice in experience good markets: Ackerberg (2003), IER
Route choice models: Fosgerau et al. (2013), Transportation Research B
Fertility and labor supply decisions: Francesconi (2002), JoLE
Residential and work-location choice: Buchinsky et al. (2014), Econometrica
Equilibrium allocations under alternative waitlist designs, deceased donor kidneys: Agarwal et al. (2021), Econometrica
Equilibrium trade in automobiles: Gillingham et al. (2022), JPE
...and many more
To this day NFXP is a very powerful method for both solving and estimating dynamic discrete choice models. In Iskhakov et al. (2016) we show that with a proper implementation Rustโs method is as powerful as the modern general-purpose solvers unleashed at the same problem today.
The components of a dynamic model โ state , decision , flow utility , motion rules, value function and policy function โ are those of Class 7, and we now fill each of them with content.
Harold Zurcherโs engine replacement problemยถ
Harold Alois Zuercher, superintendent of the Madison, Wisconsin bus company, June 16, 1926 โ June 21, 2020
States and choicesยถ
Each bus comes in for repair once a month and Zurcher chooses between ordinary maintenance and overhaul/engine replacement , so the choice set is .
The state is the mileage since the last engine overhaul. Harold observes many other attributes of the buses which come into the shop, but we focus on the main one for now.
Zurcherโs preferencesยถ
Instantaneous payoffs are given by the cost function that depends on the choice:
where is the replacement cost and the cost of maintenance with preference parameters . Replacing resets the mileage to zero, so the maintenance cost after replacement is .
Motion rules / transition probabilitiesยถ
Mileage is continuous, and the first modeling decision is how to deal with a continuous state space. Rust discretized the range of travelled miles into bins, indexed with , with . The mileage transition probability is then, for ,
Mileage in the next period can move up at most grid points, with determined from the observed distribution of mileage increments. The probabilities of an increase from any existing mileage are the same, so is a short vector.
If not replacing () the transition matrix is banded:
The last rows are adjusted so that every row sums to one: the top bin absorbs whatever mileage would fall beyond the grid. If replacing () every row is the first row of , because the mileage restarts from zero:
Dynamic optimization and the Bellman equationยถ
To minimize the discounted expected value of the costs, Zurcher should find the policy function such that maximizes
The function denotes the maximum attainable value at period ,
where is the space of policy functions , and . Using Bellmanโs principle of optimality, the value function is the solution of the functional equation
where the expectation is taken over the next period values of the state given the motion rule of the problem.
Bellman operatorยถ
The Bellman equation can be written as a fixed point equation of the Bellman operator in the functional space,
The Bellman equation is then , with the solution given by the fixed point . Everything from Class 8 on contraction mappings applies.
Classification of the Zurcher model:
Finite or infinite horizon?
Discrete or continuous choice?
Discrete or continuous states?
Which solution methods are suitable?
Making the model suitable for empirical workยถ
Zurcher observes many different attributes of the buses that come into the shop, but we as econometricians do not. Yet these are likely to be the reason for observing different behavior in the same states, that is for the same observed mileage. Without them the model predicts the same choice every time a bus shows up with a given mileage, and one observation to the contrary makes the likelihood zero.
We therefore need error terms in the model, denoted , and we are back to the random utility framework of Class 6, now in a dynamic setting.
Should the error term be part of the state space of the problem?
What would the econometrician have to observe for the model without error terms to fit the data?
Updating the Bellman equationยถ
is a new (vector) state variable, and the Bellman equation becomes
Writing the expectation out as an integral over both future states,
where is the component of the vector which corresponds to . As it stands this is a fixed point in a function of a continuous two-dimensional argument, with a double integral inside โ not something we want to evaluate thousands of times inside an estimator.
Rust assumptionsยถ
Three assumptions bring the problem back to a fixed point on the mileage grid alone.
(AS) is the same as in the static random utility models. (CI) adds to the standard independence assumptions of static RUM โ error terms independent across observations by random sampling, and independent across choices within the vector โ one assumption specific to dynamic models: conditional on , there is no serial correlation in the error terms across time. (EV) makes things simpler by providing closed form expressions for the choice probabilities. It can easily be relaxed to GEV, and with the trouble of multidimensional integration to any other continuous distribution with full support.
What the Rust assumptions allowยถ
Under (AS) and (CI) the Bellman equation reads
and the three assumptions are used in turn:
Separate out the deterministic part of the choice specific value function (AS)
Compute the expectation by parts (CI)
Use max-stability of EV1 to compute the expectation w.r.t. (EV)
The first two steps give
so that is the maximum of choice specific values plus EV1 shocks. The third step is the logsum formula of Class 6 for the expectation of that maximum:
The has been integrated out analytically, and what is left is a function of the mileage grid only.
Bellman equation in expected value function spaceยถ
Let . Then the model is the pair of equations
This is the Bellman equation in expected value function space. When the state space is discrete the integral is, of course, a simple sum over future values, and substituting into the first line gives
which defines the operator in the expected value function space,
The solution of the Bellman functional equation is a fixed point of , , and this is the fixed point we compute. It is the better one for four reasons:
is also a contraction mapping in the space of expected value functions, as shown by Ma & Stachurski (2021), so VFI is guaranteed to find the unique solution
the dimensionality of this fixed point problem is smaller than the one in value function terms, because is not included
it is numerically easier to work with the smooth expected values than with
it opens very nice numerical optimization possibilities, which is the rest of this class
Choice probabilitiesยถ
Once the fixed point is found, the optimal choice probability is given by the logit structure (assumption EV):
The choice probabilities are the basis for forming the likelihood function, and we continue with them on Thursday when talking about structural estimation.
Solution approachesยถ
The Zurcher model has the following features: infinite horizon, discretized mileage which is the only state in the formulation, so a finite state space, discrete choice, and idiosyncratic random components. Therefore the suitable solution methods are
value function iterations (VFI)
policy iterations
much faster than VFI
NewtonโKantorovich (NK) iterations = Newton-Raphson for Bellman fixed point
the reason behind the fast convergence of policy iterations
mathematically equivalent to NK iterations
Policy iterations aka Howard policy improvement algorithmยถ
Break the search of the fixed point of Bellman operator into two steps:
Policy evaluation: compute value function for fixed policy function
Policy improvement : find best policy for given value function
Iterative (back and forth) approach instead of jointly finding the value and policy functions in the fixed point problem.
Step 1: Policy evaluationยถ
For fixed policy function (decisions) solve for function :
functional equation, although simpler than Bellman equation
becomes non-linear system of equations with discrete or discretized states
becomes linear system of equations in many applications where can be expressed as matrix multiplication
Markov discrete choice problems (all randomness in transition probabilities, no interpolation)
when quadrature integration and linear interpolation are used in discretized state spaces
Step 2: Policy improvementยถ
For fixed value function find the improved optimal policy while computing:
standard evaluation of the Bellman operator
any approaches for implementation of the Bellman operator are applicable (continuous, discrete or discretized choice spaces)
Policy iterations algorithm:
1. Initialize policy function
2. Compute the value of the current policy
- by solving the Bellman equation without max operation
- assuming that current policy will be applied forever
3. Re-compute the policy function
- by applying the Bellman operator to the found value function
4. Repeat until convergence in policy and/or value function spaceWhy policy iterations? Rate of convergence!
Policy iterations are equivalent to applying Newton-Raphson method to solve for the fixed point of the Bellman operator.
NewtonโKantorovich iterationsยถ
Main idea: apply Newton-Raphson method to the fixed point equation
Kantorovich showed how to do this in functional spaces.
Write the fixed point equation as a root finding problem for the operator ,
is the identity operator
so maps a function to a difference between that function and its image under
is the zero function
The NK iteration is Newtonโs step on that equation,
is the Frรฉchet derivative of the operator
is the Frรฉchet derivative of the Bellman operator
We work with finite approximations on the discrete state space, so everything is a vector or a matrix.
Let denote the number of state points (in mileage).
is a vector of length
with the first element reused to describe the expected value of replacing
is a non-linear -valued multivariate function of
the Frรฉchet derivative is an Jacobian matrix
NK iterations on finite approximations are therefore solving a system of equations with unknowns with Newtonโs method โ see Class 5 โ in dimensions
Matrix expression for the finite approximation of the Bellman operatorยถ
On the grid, the sum over future states is a matrix product:
is an column vector
is the matrix of mileage transition probabilities
is a column vector of costs for all points in the state space, conditional on decision
is the logsum function returning an vector
the notation denotes the -th element of a vector
Only appears: replacing sends every bus to the first row, which is why the replacement value uses the single element .
Implementation of the Frรฉchet derivativeยถ
With the finite approximation of the Bellman operator
the Frรฉchet derivative w.r.t. is the matrix
so all that is needed is the derivative of the logsum function. Recall is the expectation of the maximum of with independent EV1 shocks, and let denote the corresponding choice probabilities. Differentiating with respect to any scalar ,
The derivative of the logsum is the probability-weighted derivative of its arguments. For with only the keep argument depends on it, in one row, so
where is shorthand for the probability of keeping at state point . The first column, , is special because the replacement value in every row depends on it:
where is shorthand for the probability of replacing at state point . Putting the two together, the Frรฉchet derivative in matrix notation is
The first term is with its columns scaled by the keep probabilities, the second is non-zero in the first column only, so the whole matrix is
In code this is two lines: an element-wise product of with the row vector of keep probabilities, and a correction to the first column.
NK iterations algorithmยถ
Input: model, starting point EV_0, tolerance, maxiter
Algorithm:
1. Initialize the expected value function at EV_0 (starting values matter!)
2. Apply the Bellman operator ฮ(EV_k), computing the choice probabilities
and the Frรฉchet derivative ฮ' along the way
3. NewtonโKantorovich step
EV_{k+1} = EV_k - (I - ฮ')^{-1} (I - ฮ)(EV_k)
4. Repeat from 2 until ||EV_{k+1} - EV_k|| < toleranceThe model class below follows the course architecture: the model object holds the
parameters and builds the transition matrix โ assigning n or p rebuilds the grid
and the matrix, so the object is never out of date โ bellman() is the operator with
an optional Frรฉchet derivative, and three solvers share the (maxiter, tol, callback)
signature. Every solver also accepts a starting point ev0; it defaults to zeros here,
and on Thursday the estimator will pass the previous solution instead. Parameter
values are the estimates from Rust (1987), with mileage measured in thousands and the
maintenance cost linear in it. This is the class in zurcher.py of the code
repository, which the rest of the course imports.
Source
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
class zurcher():
'''Harold Zurcher bus engine replacement model class, VFI version'''
def __init__(self,
n = 175, # number of state points
RC = 11.7257, # replacement cost
c = 2.45569, # parameter of maintance cost (theta_1)
p = [0.0937,0.4475,0.4459,0.0127], # probabilities of transitions (theta_2)
beta = 0.9999): # discount factor
'''Init for the Zurcher model object'''
self.RC, self.c, self.beta = RC, c, beta
self.p = p # transition probabilities, list; trpr is built by the n setter below
self.n = n # builds the grid and the transition matrix
@property
def n(self):
'''Attribute getter for n'''
return self.__n
@n.setter
def n(self, value):
'''Attribute n setter'''
if hasattr(self, '_zurcher__p'): # p may not be set yet when n is assigned first
assert len(self.__p) < value, 'More transition probability parameters than grid points'
self.__n = value
self.grid = np.arange(self.__n)
if hasattr(self, '_zurcher__p'):
self.trpr = self.__transition_probs()
@property
def p(self):
'''Transition probabilities theta_2, all but the last (residual) one.
Returns a copy: change p by assigning a new list, so the matrix is rebuilt'''
return list(self.__p)
@p.setter
def p(self, value):
'''Assigning p rebuilds the transition matrix'''
value = [float(v) for v in value]
assert min(value)>=0.0, 'Transition probability parameters must be non-negative'
assert sum(value)<=1.0, 'Transition probability parameters must sum up to <1'
self.__p = value
if hasattr(self, '_zurcher__n'): # n may not be set yet when p is assigned first
assert len(value) < self.__n, 'More transition probability parameters than grid points'
self.trpr = self.__transition_probs()
def __repr__(self):
'''String representation of the Zurcher model'''
return 'Rust model of bus engine replacement (id={})'.format(id(self))
def __transition_probs(self):
'''Computing the transision probability matrix'''
trpr = np.zeros((self.__n,self.__n)) # init
probs = self.p + [1-sum(self.p)] # ensure sum up to 1
for i,p in enumerate(probs):
trpr += np.diag([p]*(self.__n-i),k=i)
trpr[:,-1] = 1.-np.sum(trpr[:,:-1],axis=1) # last column absorbs the mass beyond the grid
return trpr
def bellman(self,ev0,deriv=False):
'''Bellman operator for the model
Depending on deriv argument, returns 2 or 3 outputs (Frรฉchet derivative)
'''
x = self.grid # points in the next period state
mcost = -0.001*x*self.c # 1-dim array of maintenance costs
vx0 = mcost + self.beta * ev0 # 1-dim array v(x,0), keep
vx1 = mcost[0] - self.RC + self.beta * ev0[0] # 1-dim array v(x,1), replace
M = np.maximum(vx0,vx1) # de-max values to avoid exp(large number)
logsum = M + np.log(np.exp(vx0-M) + np.exp(vx1-M))
ev1 = self.trpr @ logsum # 1-dim array after matrix multiplication
pk = 1/( np.exp(vx1-vx0)+1 ) # choice prob to keep
if not deriv:
return ev1, pk
# Frรฉchet derivative
dev1 = self.beta * self.trpr * pk[np.newaxis,:] # element-wise, pk in rows
dev1[:,0] += self.beta * self.trpr @ (1-pk) # w.r.t. EV[0] special case
return ev1, pk, dev1
def _start(self, ev0):
'''Starting point for the solvers: zeros unless one is given'''
return np.zeros(self.n) if ev0 is None else np.asarray(ev0, dtype=float).copy()
def solve_vfi(self, maxiter=100, tol=1e-6, callback=None, ev0=None):
'''Solves the Rust model using value function iterations
'''
ev0 = self._start(ev0) # initial point for VFI
err0 = 1.0 # initial lagged error
for iter in range(maxiter): # main loop
ev1, pk = self.bellman(ev0) # update approximation
err = np.amax(np.abs(ev0-ev1))
if callback:
callback(iter=iter,model=self,ev1=ev1,ev0=ev0,err=err,err_prev=err0,pk=pk,method='vfi',itertype='sa')
if err<tol:
break # break out if converged
ev0 = ev1 # get ready to the next iteration
err0 = err
else:
raise RuntimeError('Failed to converge in %d iterations'%maxiter)
return ev1, pk
def solve_nk(self, maxiter=100, tol=1e-6, callback=None, ev0=None):
'''Solves the model using the Newton-Kantorovich iterations
'''
ev0 = self._start(ev0) # initial point
err0 = 1.0 # initial lagged error
for iter in range(maxiter):
ev1,pk,dev = self.bellman(ev0,deriv=True) # compute with Frรฉchet derivative
ev1 = ev0 - np.linalg.solve(np.eye(self.n)-dev,ev0 - ev1) # NK step
err = np.max(np.abs(ev1-ev0))
if callback:
callback(iter=iter,model=self,ev1=ev1,ev0=ev0,err=err,err_prev=err0,pk=pk,method='nk',itertype='nk')
if err < tol:
break # break out if converged
ev0 = ev1 # get ready to the next iteration
err0 = err
else:
raise RuntimeError('Failed to converge in %d iterations'%maxiter)
ev1,pk = self.bellman(ev1) # compute choice probabilities after convergence
return ev1,pk
def solve_poly(self,
maxiter=100,
tol=1e-10,
sa_min=5, # minimum number of contraction steps
sa_max=25, # maximum number of contraction steps
switch_tol=0.025, # tolerance of the switching rule
callback=None,
ev0=None):
'''Solves the model using the poly-algorithm'''
ev0 = self._start(ev0) # initial point
err0 = 1.0 # initial lagged error
nk = False # start with successive approximations
for iter in range(maxiter):
ev1,pk,dev = self.bellman(ev0,deriv=True) # update EV for both types of a step
err = np.max(np.abs(ev1-ev0))
nk = True if iter>= sa_max else nk # have to switch to NK after sa_max
nk = nk or (iter>=sa_min and abs(err/err0 - self.beta)<switch_tol) # check if need to switch to NK
if nk:
ev1 = ev0 - np.linalg.solve(np.eye(self.n)-dev,ev0 - ev1) # NK step
err = np.max(np.abs(ev1-ev0))
if callback:
itertype = 'nk' if nk else 'sa' # label for the iteration type
callback(iter=iter,model=self,ev1=ev1,ev0=ev0,err=err,err_prev=err0,pk=pk,method='poly',itertype=itertype)
if err < tol:
break # break out if converged
ev0 = ev1 # get ready to the next iteration
err0 = err
else:
raise RuntimeError('No convergence: maximum number of iterations achieved! Increase maxiter')
ev1,pk = self.bellman(ev1) # compute choice probabilities after convergence
return ev1,pk
def solve_show(self,solver='vfi',verbosity=0,plot=True,**kvargs):
'''Illustrate solution for given solver = {vfi,nk,poly} and
print errors/relative errors from iterations (when verbose=True)
All other arguments are passed to the solver
'''
if solver=='vfi':
chosen_solver = self.solve_vfi
elif solver=='nk':
chosen_solver = self.solve_nk
elif solver=='poly':
chosen_solver = self.solve_poly
else:
raise RuntimeError('Unknown solver in solve_show()')
if plot:
fig1, (ax1,ax2) = plt.subplots(1,2,figsize=(14,8))
ax1.grid(visible=True, which='both', color='0.65', linestyle='-')
ax2.grid(visible=True, which='both', color='0.65', linestyle='-')
ax1.set_xlabel('Mileage grid')
ax2.set_xlabel('Mileage grid')
ax1.set_title(f'Value function ({solver})')
ax2.set_title(f'Probability of replacing the engine ({solver})')
def callback(**argvars):
iter,itertype,err,derr = argvars['iter'],argvars['itertype'],argvars['err'],argvars['err_prev']
mod, ev, pk = argvars['model'],argvars['ev1'],argvars['pk']
if verbosity>1:
if iter==0:
print('Solver = %s'%solver)
print('-'*42)
print('%7s %16s %16s'%('iter','err','err(i)/err(i-1)'))
print('-'*42)
print('%4d %2s %16.4e %16.12f'%(iter,itertype[:2],err,err/derr))
elif verbosity>0:
if iter==0:
print('Solver = %s'%solver)
print('-'*22)
print('%4s %16s'%('iter','err'))
print('-'*22)
print('%4d %16.4e'%(iter,err))
if plot:
ax1.plot(mod.grid,ev,color='k',alpha=0.25)
ax2.plot(mod.grid,pk,color='k',alpha=0.25)
callback.nriter = iter+1 # number of iterations run, saved in function object attribute
# run the chosen solver
ev,pk = chosen_solver(callback=callback,**kvargs)
if plot:
# add solutions
ax1.plot(self.grid,ev,color='r',linewidth=2.5)
ax2.plot(self.grid,pk,color='r',linewidth=2.5)
plt.show()
print('{} solved with {} in {} iterations'.format(self,solver,callback.nriter))
return ev,pkThe two solvers find the same fixed point, and the grey lines show how differently they get there:
# compare SA, NK
model = zurcher(beta=0.9) # try different value of beta
ev1,pk1 = model.solve_show(maxiter=1500)
ev2,pk2 = model.solve_show(solver='nk')
print()
print('Max diff between value functions is ' ,np.amax(np.abs(ev1-ev2)))
print('Max diff between policy functions is',np.amax(np.abs(pk1-pk2)))
Rust model of bus engine replacement (id=140628910063152) solved with vfi in 124 iterations
Rust model of bus engine replacement (id=140628910063152) solved with nk in 3 iterations
Max diff between value functions is 8.95709148851509e-06
Max diff between policy functions is 9.113154675333135e-12
Does the VFI algorithm always converge?
What determines the speed of convergence of the VFI algorithm?
Does the NK algorithm always converge?
Properties of VFI vs NewtonโKantorovichยถ
VFI is globally convergent, because the Bellman operator is a contraction mapping with a single fixed point, but its convergence rate is : very slow in approaching the fixed point when is close to one. NewtonโKantorovich has quadratic convergence but is sensitive to the starting point, like every Newton method. Raising from 0.9 to 0.975 makes the difference visible in the iteration counts:
# compare SA, NK
model = zurcher(beta=0.975)
ev1,pk1 = model.solve_show(maxiter=1500,verbosity=0,plot=False)
ev2,pk2 = model.solve_show(solver='nk',verbosity=0,plot=False)
print()
print('Max diff between value functions is ' ,np.amax(np.abs(ev1-ev2)))
print('Max diff between policy functions is',np.amax(np.abs(pk1-pk2)))Rust model of bus engine replacement (id=140627787649616) solved with vfi in 497 iterations
Rust model of bus engine replacement (id=140627787649616) solved with nk in 8 iterations
Max diff between value functions is 3.884109376706135e-05
Max diff between policy functions is 2.4160659028993337e-09
Poly-algorithmยถ
The NK method may not be convergent at the initial point, whereas successive approximation (SA) iterations are always convergent. The poly-algorithm is the combination of the two: start with SA iterations, and at approximately the optimal time switch to NK iterations. The question is what โoptimal timeโ means.
When to switch to NK iterations?ยถ
Suppose the current approximation is a constant away from the fixed point, . Then two consecutive SA errors are
and the ratio of the two errors is exactly when the current approximation is a constant away from the fixed point. An NK iteration โstrips awayโ such a constant in one step. Thus, switch to NK iterations when is close to . The ratio is printed in the third column below:
# when to switch from SA to NK
model = zurcher(beta=0.975)
model.solve_show(maxiter=1500,verbosity=2,plot=False);Solver = vfi
------------------------------------------
iter err err(i)/err(i-1)
------------------------------------------
0 sa 4.2728e-01 0.427277667463
1 sa 4.1659e-01 0.974985156037
2 sa 4.0617e-01 0.974977898685
3 sa 3.9600e-01 0.974967530759
4 sa 3.8608e-01 0.974952914413
5 sa 3.7640e-01 0.974932594839
6 sa 3.6696e-01 0.974905009462
7 sa 3.5773e-01 0.974867836074
8 sa 3.4873e-01 0.974819168371
9 sa 3.3992e-01 0.974756037945
10 sa 3.3131e-01 0.974674240553
11 sa 3.2289e-01 0.974575376329
12 sa 3.1464e-01 0.974446349105
13 sa 3.0655e-01 0.974300520436
14 sa 2.9862e-01 0.974112724869
15 sa 2.9082e-01 0.973904920662
16 sa 2.8316e-01 0.973656096106
17 sa 2.7562e-01 0.973365083189
18 sa 2.6820e-01 0.973069628117
19 sa 2.6088e-01 0.972704596040
20 sa 2.5366e-01 0.972323312562
21 sa 2.4654e-01 0.971949651340
22 sa 2.3952e-01 0.971507821865
23 sa 2.3258e-01 0.971044307238
24 sa 2.2576e-01 0.970670088415
25 sa 2.1904e-01 0.970230588772
26 sa 2.1242e-01 0.969791221048
27 sa 2.0592e-01 0.969373726567
28 sa 1.9954e-01 0.969042250559
29 sa 1.9331e-01 0.968766942633
30 sa 1.8722e-01 0.968482332963
31 sa 1.8127e-01 0.968253922437
32 sa 1.7549e-01 0.968087699003
33 sa 1.6987e-01 0.967987211283
34 sa 1.6443e-01 0.967953475795
35 sa 1.5916e-01 0.967985087030
36 sa 1.5408e-01 0.968078487981
37 sa 1.4919e-01 0.968228349290
38 sa 1.4448e-01 0.968428005047
39 sa 1.3995e-01 0.968669899577
40 sa 1.3560e-01 0.968946009560
41 sa 1.3143e-01 0.969248216930
42 sa 1.2743e-01 0.969568618533
43 sa 1.2360e-01 0.969899767409
44 sa 1.1992e-01 0.970234847204
45 sa 1.1639e-01 0.970567785770
46 sa 1.1300e-01 0.970893316552
47 sa 1.0975e-01 0.971206997381
48 sa 1.0662e-01 0.971505196210
49 sa 1.0361e-01 0.971785052539
50 sa 1.0072e-01 0.972044422082
51 sa 9.7925e-02 0.972281810899
52 sa 9.5235e-02 0.972524397195
53 sa 9.2641e-02 0.972764776593
54 sa 9.0134e-02 0.972944023297
55 sa 8.7710e-02 0.973099593423
56 sa 8.5366e-02 0.973277772047
57 sa 8.3101e-02 0.973464140669
58 sa 8.0904e-02 0.973561203087
59 sa 7.8777e-02 0.973712319965
60 sa 7.6716e-02 0.973842234165
61 sa 7.4714e-02 0.973898288277
62 sa 7.2777e-02 0.974076641241
63 sa 7.0891e-02 0.974091073911
64 sa 6.9065e-02 0.974242653413
65 sa 6.7287e-02 0.974255178775
66 sa 6.5563e-02 0.974379430258
67 sa 6.3884e-02 0.974382651720
68 sa 6.2255e-02 0.974499936526
69 sa 6.0666e-02 0.974477863212
70 sa 5.9125e-02 0.974605561245
71 sa 5.7622e-02 0.974567831503
72 sa 5.6162e-02 0.974675310244
73 sa 5.4739e-02 0.974656991615
74 sa 5.3355e-02 0.974710686210
75 sa 5.2007e-02 0.974735044237
76 sa 5.0692e-02 0.974727589394
77 sa 4.9415e-02 0.974803664879
78 sa 4.8168e-02 0.974755914814
79 sa 4.6956e-02 0.974837929159
80 sa 4.5773e-02 0.974814952630
81 sa 4.4621e-02 0.974822906957
82 sa 4.3499e-02 0.974867882447
83 sa 4.2404e-02 0.974815613596
84 sa 4.1339e-02 0.974899405624
85 sa 4.0300e-02 0.974862755710
86 sa 3.9287e-02 0.974865780117
87 sa 3.8302e-02 0.974906196811
88 sa 3.7338e-02 0.974851231802
89 sa 3.6402e-02 0.974923045769
90 sa 3.5488e-02 0.974891499082
91 sa 3.4596e-02 0.974879387597
92 sa 3.3729e-02 0.974929612245
93 sa 3.2882e-02 0.974873189754
94 sa 3.2057e-02 0.974926474993
95 sa 3.1253e-02 0.974909746507
96 sa 3.0468e-02 0.974877825399
97 sa 2.9704e-02 0.974945056120
98 sa 2.8958e-02 0.974887934600
99 sa 2.8232e-02 0.974919890745
100 sa 2.7524e-02 0.974922606917
101 sa 2.6832e-02 0.974868989623
102 sa 2.6160e-02 0.974956523865
103 sa 2.5504e-02 0.974899132392
104 sa 2.4864e-02 0.974908750625
105 sa 2.4241e-02 0.974932896779
106 sa 2.3632e-02 0.974874944979
107 sa 2.3040e-02 0.974948184810
108 sa 2.2461e-02 0.974908720068
109 sa 2.1898e-02 0.974895759242
110 sa 2.1349e-02 0.974942061979
111 sa 2.0813e-02 0.974884147137
112 sa 2.0291e-02 0.974934615992
113 sa 1.9782e-02 0.974917651285
114 sa 1.9276e-02 0.974442324579
115 sa 1.8768e-02 0.973647615683
116 sa 1.8258e-02 0.972800995844
117 sa 1.7745e-02 0.971899500846
118 sa 1.7229e-02 0.970941912168
119 sa 1.6711e-02 0.969929803275
120 sa 1.6191e-02 0.968868531423
121 sa 1.5669e-02 0.967767969539
122 sa 1.5172e-02 0.968243679246
123 sa 1.4790e-02 0.974879134735
124 sa 1.4409e-02 0.974190667298
125 sa 1.4027e-02 0.973536217403
126 sa 1.3648e-02 0.972954555652
127 sa 1.3272e-02 0.972469171131
128 sa 1.2902e-02 0.972092733821
129 sa 1.2538e-02 0.971814669992
130 sa 1.2183e-02 0.971631827049
131 sa 1.1836e-02 0.971530682589
132 sa 1.1498e-02 0.971501938002
133 sa 1.1171e-02 0.971523275621
134 sa 1.0854e-02 0.971593960300
135 sa 1.0546e-02 0.971696642005
136 sa 1.0249e-02 0.971822242956
137 sa 9.9619e-03 0.971964481885
138 sa 9.6842e-03 0.972123123013
139 sa 9.4159e-03 0.972288457325
140 sa 9.1565e-03 0.972451170140
141 sa 8.9057e-03 0.972614931055
142 sa 8.6634e-03 0.972787714695
143 sa 8.4290e-03 0.972942421442
144 sa 8.2021e-03 0.973087799203
145 sa 7.9827e-03 0.973248816148
146 sa 7.7703e-03 0.973388631869
147 sa 7.5645e-03 0.973518934386
148 sa 7.3651e-03 0.973642278530
149 sa 7.1718e-03 0.973757211091
150 sa 6.9844e-03 0.973862620201
151 sa 6.8026e-03 0.973975204398
152 sa 6.6262e-03 0.974072553481
153 sa 6.4550e-03 0.974152310416
154 sa 6.2886e-03 0.974220630047
155 sa 6.1270e-03 0.974313921035
156 sa 5.9700e-03 0.974376033104
157 sa 5.8173e-03 0.974418139010
158 sa 5.6690e-03 0.974507447965
159 sa 5.5246e-03 0.974532451652
160 sa 5.3843e-03 0.974594759405
161 sa 5.2477e-03 0.974625994713
162 sa 5.1147e-03 0.974670050757
163 sa 4.9853e-03 0.974699222394
164 sa 4.8594e-03 0.974735343938
165 sa 4.7367e-03 0.974754402472
166 sa 4.6173e-03 0.974791571133
167 sa 4.5009e-03 0.974793950852
168 sa 4.3877e-03 0.974839827553
169 sa 4.2772e-03 0.974827574342
170 sa 4.1697e-03 0.974873930177
171 sa 4.0650e-03 0.974866852219
172 sa 3.9629e-03 0.974885637004
173 sa 3.8634e-03 0.974900748082
174 sa 3.7664e-03 0.974889016101
175 sa 3.6720e-03 0.974930246334
176 sa 3.5798e-03 0.974912175141
177 sa 3.4901e-03 0.974929779093
178 sa 3.4026e-03 0.974937031454
179 sa 3.3173e-03 0.974919583096
180 sa 3.2342e-03 0.974959317818
181 sa 3.1532e-03 0.974939050705
182 sa 3.0742e-03 0.974946267103
183 sa 2.9972e-03 0.974958716208
184 sa 2.9221e-03 0.974937686756
185 sa 2.8489e-03 0.974968002465
186 sa 2.7776e-03 0.974955531852
187 sa 2.7080e-03 0.974947950285
188 sa 2.6402e-03 0.974972440676
189 sa 2.5741e-03 0.974950600776
190 sa 2.5096e-03 0.974964742321
191 sa 2.4468e-03 0.974966597118
192 sa 2.3855e-03 0.974944523624
193 sa 2.3258e-03 0.974979946654
194 sa 2.2676e-03 0.974959927591
195 sa 2.2108e-03 0.974956824710
196 sa 2.1555e-03 0.974975050782
197 sa 2.1015e-03 0.974952748669
198 sa 2.0489e-03 0.974970626299
199 sa 1.9976e-03 0.974967607221
200 sa 1.9475e-03 0.974946933934
201 sa 1.8988e-03 0.974982348120
202 sa 1.8513e-03 0.974959968761
203 sa 1.8049e-03 0.974960184441
204 sa 1.7598e-03 0.974974599124
205 sa 1.7157e-03 0.974952232589
206 sa 1.6727e-03 0.974973271840
207 sa 1.6309e-03 0.974966797135
208 sa 1.5900e-03 0.974944457876
209 sa 1.5498e-03 0.974701813459
210 sa 1.5101e-03 0.974407167329
211 sa 1.4710e-03 0.974115684452
212 sa 1.4325e-03 0.973828205861
213 sa 1.3946e-03 0.973545598477
214 sa 1.3573e-03 0.973268752638
215 sa 1.3207e-03 0.972998579032
216 sa 1.2847e-03 0.972736005191
217 sa 1.2493e-03 0.972481971466
218 sa 1.2147e-03 0.972237426465
219 sa 1.1814e-03 0.972637892569
220 sa 1.1518e-03 0.974951217352
221 sa 1.1227e-03 0.974712685543
222 sa 1.0941e-03 0.974495690995
223 sa 1.0660e-03 0.974309379414
224 sa 1.0384e-03 0.974155843925
225 sa 1.0114e-03 0.974034584402
226 sa 9.8509e-04 0.973943128620
227 sa 9.5936e-04 0.973878016459
228 sa 9.3426e-04 0.973835249673
229 sa 9.0979e-04 0.973810313460
230 sa 8.8595e-04 0.973799322791
231 sa 8.6275e-04 0.973806306198
232 sa 8.4016e-04 0.973819901434
233 sa 8.1818e-04 0.973842632946
234 sa 7.9681e-04 0.973873124117
235 sa 7.7602e-04 0.973908854307
236 sa 7.5580e-04 0.973947234364
237 sa 7.3614e-04 0.973991048395
238 sa 7.1703e-04 0.974038393599
239 sa 6.9845e-04 0.974081910804
240 sa 6.8038e-04 0.974130880202
241 sa 6.6281e-04 0.974180726386
242 sa 6.4573e-04 0.974225575142
243 sa 6.2911e-04 0.974271409586
244 sa 6.1296e-04 0.974317441612
245 sa 5.9724e-04 0.974362939950
246 sa 5.8196e-04 0.974407238039
247 sa 5.6709e-04 0.974449740273
248 sa 5.5262e-04 0.974489926507
249 sa 5.3854e-04 0.974527354758
250 sa 5.2485e-04 0.974561662037
251 sa 5.1151e-04 0.974592563528
252 sa 4.9853e-04 0.974631024089
253 sa 4.8590e-04 0.974662289687
254 sa 4.7360e-04 0.974684544694
255 sa 4.6163e-04 0.974714660413
256 sa 4.4997e-04 0.974742145135
257 sa 4.3861e-04 0.974754985440
258 sa 4.2755e-04 0.974790768691
259 sa 4.1678e-04 0.974799416103
260 sa 4.0628e-04 0.974824563740
261 sa 3.9606e-04 0.974836188537
262 sa 3.8610e-04 0.974855080567
263 sa 3.7640e-04 0.974865638009
264 sa 3.6694e-04 0.974882285106
265 sa 3.5773e-04 0.974888311176
266 sa 3.4875e-04 0.974906305202
267 sa 3.4000e-04 0.974904887595
268 sa 3.3147e-04 0.974927379527
269 sa 3.2316e-04 0.974921704812
270 sa 3.1506e-04 0.974940221539
271 sa 3.0717e-04 0.974939291038
272 sa 2.9947e-04 0.974945397979
273 sa 2.9197e-04 0.974954676860
274 sa 2.8466e-04 0.974947249734
275 sa 2.7753e-04 0.974967906728
276 sa 2.7058e-04 0.974960170066
277 sa 2.6381e-04 0.974965566489
278 sa 2.5720e-04 0.974971658481
279 sa 2.5077e-04 0.974963024558
280 sa 2.4449e-04 0.974979998305
281 sa 2.3837e-04 0.974972972415
282 sa 2.3241e-04 0.974973210967
283 sa 2.2659e-04 0.974982088032
284 sa 2.2092e-04 0.974972721684
285 sa 2.1539e-04 0.974983054548
286 sa 2.1001e-04 0.974980967230
287 sa 2.0475e-04 0.974973885721
288 sa 1.9963e-04 0.974988777168
289 sa 1.9463e-04 0.974979011792
290 sa 1.8976e-04 0.974981363660
291 sa 1.8502e-04 0.974986371948
292 sa 1.8039e-04 0.974976503352
293 sa 1.7588e-04 0.974988066189
294 sa 1.7148e-04 0.974983567317
295 sa 1.6719e-04 0.974977415252
296 sa 1.6300e-04 0.974990508552
297 sa 1.5893e-04 0.974980516945
298 sa 1.5495e-04 0.974983445115
299 sa 1.5107e-04 0.974987321892
300 sa 1.4729e-04 0.974977324114
301 sa 1.4361e-04 0.974989268680
302 sa 1.4002e-04 0.974984043729
303 sa 1.3651e-04 0.974978281751
304 sa 1.3310e-04 0.974990744919
305 sa 1.2977e-04 0.974980721653
306 sa 1.2652e-04 0.974983971571
307 sa 1.2336e-04 0.974987383899
308 sa 1.2027e-04 0.974977384312
309 sa 1.1725e-04 0.974875056333
310 sa 1.1429e-04 0.974744215194
311 sa 1.1139e-04 0.974616422921
312 sa 1.0855e-04 0.974492048880
313 sa 1.0576e-04 0.974371450713
314 sa 1.0304e-04 0.974254973773
315 sa 1.0038e-04 0.974142950394
316 sa 9.7771e-05 0.974035699006
317 sa 9.5223e-05 0.973933523453
318 sa 9.2731e-05 0.973836711944
319 sa 9.0297e-05 0.973745536287
320 sa 8.8017e-05 0.974748877580
321 sa 8.5809e-05 0.974916964960
322 sa 8.3648e-05 0.974820564926
323 sa 8.1535e-05 0.974736522706
324 sa 7.9470e-05 0.974666699492
325 sa 7.7452e-05 0.974611257558
326 sa 7.5482e-05 0.974568879462
327 sa 7.3560e-05 0.974537709974
328 sa 7.1686e-05 0.974517191892
329 sa 6.9858e-05 0.974504438039
330 sa 6.8077e-05 0.974497754111
331 sa 6.6341e-05 0.974499171361
332 sa 6.4649e-05 0.974503426588
333 sa 6.3001e-05 0.974511835365
334 sa 6.1396e-05 0.974523669667
335 sa 5.9833e-05 0.974537878310
336 sa 5.8310e-05 0.974553390278
337 sa 5.6828e-05 0.974570639189
338 sa 5.5384e-05 0.974590704367
339 sa 5.3977e-05 0.974607695595
340 sa 5.2608e-05 0.974629235178
341 sa 5.1274e-05 0.974648523239
342 sa 4.9975e-05 0.974667072960
343 sa 4.8710e-05 0.974686122398
344 sa 4.7478e-05 0.974706543515
345 sa 4.6278e-05 0.974725349592
346 sa 4.5109e-05 0.974742998273
347 sa 4.3971e-05 0.974759965876
348 sa 4.2862e-05 0.974777732717
349 sa 4.1781e-05 0.974794547332
350 sa 4.0729e-05 0.974809469470
351 sa 3.9703e-05 0.974823048655
352 sa 3.8704e-05 0.974835186022
353 sa 3.7731e-05 0.974850569319
354 sa 3.6783e-05 0.974863240946
355 sa 3.5858e-05 0.974871776666
356 sa 3.4958e-05 0.974884718477
357 sa 3.4080e-05 0.974894807300
358 sa 3.3225e-05 0.974900520883
359 sa 3.2391e-05 0.974914706606
360 sa 3.1579e-05 0.974917201567
361 sa 3.0787e-05 0.974929390832
362 sa 3.0015e-05 0.974931820436
363 sa 2.9263e-05 0.974941752257
364 sa 2.8530e-05 0.974943509515
365 sa 2.7815e-05 0.974952799030
366 sa 2.7119e-05 0.974952490870
367 sa 2.6440e-05 0.974962584313
368 sa 2.5778e-05 0.974960488398
369 sa 2.5132e-05 0.974969749005
370 sa 2.4503e-05 0.974968736428
371 sa 2.3890e-05 0.974973471621
372 sa 2.3292e-05 0.974975970608
373 sa 2.2709e-05 0.974975446428
374 sa 2.2141e-05 0.974982326499
375 sa 2.1587e-05 0.974979168765
376 sa 2.1047e-05 0.974984727672
377 sa 2.0521e-05 0.974984536199
378 sa 2.0007e-05 0.974983663480
379 sa 1.9507e-05 0.974989325384
380 sa 1.9019e-05 0.974985685075
381 sa 1.8543e-05 0.974989664757
382 sa 1.8080e-05 0.974989848382
383 sa 1.7627e-05 0.974986783621
384 sa 1.7187e-05 0.974993671496
385 sa 1.6757e-05 0.974989733800
386 sa 1.6338e-05 0.974990900750
387 sa 1.5929e-05 0.974993198590
388 sa 1.5531e-05 0.974989166240
389 sa 1.5142e-05 0.974994353232
390 sa 1.4764e-05 0.974992381581
391 sa 1.4394e-05 0.974990171512
392 sa 1.4034e-05 0.974995478774
393 sa 1.3684e-05 0.974991335108
394 sa 1.3341e-05 0.974992981781
395 sa 1.3008e-05 0.974994307793
396 sa 1.2682e-05 0.974990142086
397 sa 1.2365e-05 0.974995586378
398 sa 1.2056e-05 0.974993033990
399 sa 1.1755e-05 0.974991026720
400 sa 1.1461e-05 0.974995897405
401 sa 1.1174e-05 0.974991700915
402 sa 1.0895e-05 0.974993466244
403 sa 1.0622e-05 0.974994527494
404 sa 1.0357e-05 0.974990335950
405 sa 1.0098e-05 0.974995853916
406 sa 9.8451e-06 0.974993139368
407 sa 9.5988e-06 0.974991226231
408 sa 9.3588e-06 0.974995941323
409 sa 9.1248e-06 0.974991745058
410 sa 8.8962e-06 0.974952086966
411 sa 8.6729e-06 0.974897239640
412 sa 8.4547e-06 0.974843773289
413 sa 8.2416e-06 0.974791859529
414 sa 8.0334e-06 0.974741662671
415 sa 7.8301e-06 0.974693336655
416 sa 7.6316e-06 0.974647028288
417 sa 7.4378e-06 0.974602874501
418 sa 7.2486e-06 0.974561002135
419 sa 7.0639e-06 0.974521528350
420 sa 6.8837e-06 0.974484559030
421 sa 6.7106e-06 0.974862258232
422 sa 6.5426e-06 0.974968541943
423 sa 6.3786e-06 0.974928565083
424 sa 6.2185e-06 0.974893722550
425 sa 6.0622e-06 0.974864727264
426 sa 5.9097e-06 0.974841612212
427 sa 5.7609e-06 0.974823964495
428 sa 5.6158e-06 0.974811151994
429 sa 5.4743e-06 0.974802430813
430 sa 5.3363e-06 0.974796974437
431 sa 5.2018e-06 0.974794871140
432 sa 5.0707e-06 0.974794489633
433 sa 4.9429e-06 0.974797032923
434 sa 4.8183e-06 0.974800415726
435 sa 4.6969e-06 0.974805179163
436 sa 4.5786e-06 0.974810892368
437 sa 4.4633e-06 0.974817263801
438 sa 4.3509e-06 0.974825358145
439 sa 4.2414e-06 0.974832167941
440 sa 4.1347e-06 0.974840568667
441 sa 4.0307e-06 0.974847641207
442 sa 3.9294e-06 0.974857102701
443 sa 3.8306e-06 0.974864618519
444 sa 3.7344e-06 0.974872323502
445 sa 3.6406e-06 0.974880091019
446 sa 3.5491e-06 0.974887798074
447 sa 3.4600e-06 0.974895334378
448 sa 3.3732e-06 0.974902603179
449 sa 3.2886e-06 0.974909513573
450 sa 3.2061e-06 0.974915991701
451 sa 3.1257e-06 0.974921976102
452 sa 3.0473e-06 0.974927413659
453 sa 2.9709e-06 0.974933544620
454 sa 2.8965e-06 0.974939500305
455 sa 2.8239e-06 0.974943545984
456 sa 2.7531e-06 0.974947788273
457 sa 2.6842e-06 0.974953690333
458 sa 2.6170e-06 0.974956190516
459 sa 2.5514e-06 0.974960890370
460 sa 2.4876e-06 0.974964174717
461 sa 2.4253e-06 0.974966771118
462 sa 2.3646e-06 0.974970885403
463 sa 2.3054e-06 0.974972126985
464 sa 2.2477e-06 0.974976358794
465 sa 2.1915e-06 0.974976940326
466 sa 2.1366e-06 0.974980672088
467 sa 2.0832e-06 0.974981221730
468 sa 2.0311e-06 0.974983925976
469 sa 1.9803e-06 0.974985000440
470 sa 1.9307e-06 0.974986237696
471 sa 1.8824e-06 0.974988321579
472 sa 1.8354e-06 0.974987724906
473 sa 1.7895e-06 0.974991237037
474 sa 1.7447e-06 0.974990065522
475 sa 1.7011e-06 0.974992246161
476 sa 1.6585e-06 0.974992515500
477 sa 1.6171e-06 0.974992262401
478 sa 1.5766e-06 0.974994683642
479 sa 1.5372e-06 0.974993280366
480 sa 1.4988e-06 0.974995164061
481 sa 1.4613e-06 0.974995140248
482 sa 1.4247e-06 0.974994259016
483 sa 1.3891e-06 0.974996824577
484 sa 1.3544e-06 0.974995273525
485 sa 1.3205e-06 0.974996192082
486 sa 1.2875e-06 0.974996774726
487 sa 1.2553e-06 0.974995173631
488 sa 1.2239e-06 0.974997768037
489 sa 1.1933e-06 0.974996541837
490 sa 1.1635e-06 0.974996168444
491 sa 1.1344e-06 0.974997845325
492 sa 1.1060e-06 0.974996184116
493 sa 1.0784e-06 0.974997393234
494 sa 1.0514e-06 0.974997419250
495 sa 1.0251e-06 0.974995744736
496 sa 9.9949e-07 0.974998498597
Rust model of bus engine replacement (id=140628238212832) solved with vfi in 497 iterations
With a tight switching tolerance the poly-algorithm spends a few SA steps, then a few NK steps, and beats both pure solvers:
# compare SA, NK and polyalgorithm
m = zurcher(beta=0.975)
ev,pk = m.solve_show(tol=1e-10,maxiter=1500)
ev,pk = m.solve_show(tol=1e-10,solver='nk',plot=False)
polyset = {'sa_min':10,
'sa_max':100,
'switch_tol':0.000215,
}
ev,pk = m.solve_show(tol=1e-10,verbosity=2,solver='poly',**polyset)
Rust model of bus engine replacement (id=140628910063488) solved with vfi in 861 iterations
Rust model of bus engine replacement (id=140628910063488) solved with nk in 8 iterations
Solver = poly
------------------------------------------
iter err err(i)/err(i-1)
------------------------------------------
0 sa 4.2728e-01 0.427277667463
1 sa 4.1659e-01 0.974985156037
2 sa 4.0617e-01 0.974977898685
3 sa 3.9600e-01 0.974967530759
4 sa 3.8608e-01 0.974952914413
5 sa 3.7640e-01 0.974932594839
6 sa 3.6696e-01 0.974905009462
7 sa 3.5773e-01 0.974867836074
8 sa 3.4873e-01 0.974819168371
9 sa 3.3992e-01 0.974756037945
10 sa 3.3131e-01 0.974674240553
11 sa 3.2289e-01 0.974575376329
12 sa 3.1464e-01 0.974446349105
13 sa 3.0655e-01 0.974300520436
14 sa 2.9862e-01 0.974112724869
15 sa 2.9082e-01 0.973904920662
16 sa 2.8316e-01 0.973656096106
17 sa 2.7562e-01 0.973365083189
18 sa 2.6820e-01 0.973069628117
19 sa 2.6088e-01 0.972704596040
20 sa 2.5366e-01 0.972323312562
21 sa 2.4654e-01 0.971949651340
22 sa 2.3952e-01 0.971507821865
23 sa 2.3258e-01 0.971044307238
24 sa 2.2576e-01 0.970670088415
25 sa 2.1904e-01 0.970230588772
26 sa 2.1242e-01 0.969791221048
27 sa 2.0592e-01 0.969373726567
28 sa 1.9954e-01 0.969042250559
29 sa 1.9331e-01 0.968766942633
30 sa 1.8722e-01 0.968482332963
31 sa 1.8127e-01 0.968253922437
32 sa 1.7549e-01 0.968087699003
33 sa 1.6987e-01 0.967987211283
34 sa 1.6443e-01 0.967953475795
35 sa 1.5916e-01 0.967985087030
36 sa 1.5408e-01 0.968078487981
37 sa 1.4919e-01 0.968228349290
38 sa 1.4448e-01 0.968428005047
39 sa 1.3995e-01 0.968669899577
40 sa 1.3560e-01 0.968946009560
41 sa 1.3143e-01 0.969248216930
42 sa 1.2743e-01 0.969568618533
43 sa 1.2360e-01 0.969899767409
44 sa 1.1992e-01 0.970234847204
45 sa 1.1639e-01 0.970567785770
46 sa 1.1300e-01 0.970893316552
47 sa 1.0975e-01 0.971206997381
48 sa 1.0662e-01 0.971505196210
49 sa 1.0361e-01 0.971785052539
50 sa 1.0072e-01 0.972044422082
51 sa 9.7925e-02 0.972281810899
52 sa 9.5235e-02 0.972524397195
53 sa 9.2641e-02 0.972764776593
54 sa 9.0134e-02 0.972944023297
55 sa 8.7710e-02 0.973099593423
56 sa 8.5366e-02 0.973277772047
57 sa 8.3101e-02 0.973464140669
58 sa 8.0904e-02 0.973561203087
59 sa 7.8777e-02 0.973712319965
60 sa 7.6716e-02 0.973842234165
61 sa 7.4714e-02 0.973898288277
62 sa 7.2777e-02 0.974076641241
63 sa 7.0891e-02 0.974091073911
64 sa 6.9065e-02 0.974242653413
65 sa 6.7287e-02 0.974255178775
66 sa 6.5563e-02 0.974379430258
67 sa 6.3884e-02 0.974382651720
68 sa 6.2255e-02 0.974499936526
69 sa 6.0666e-02 0.974477863212
70 sa 5.9125e-02 0.974605561245
71 sa 5.7622e-02 0.974567831503
72 sa 5.6162e-02 0.974675310244
73 sa 5.4739e-02 0.974656991615
74 sa 5.3355e-02 0.974710686210
75 sa 5.2007e-02 0.974735044237
76 sa 5.0692e-02 0.974727589394
77 nk 1.7465e+00 34.452068933883
78 nk 6.4457e-03 0.003690746725
79 nk 3.5956e-06 0.000557818120
80 nk 1.7195e-12 0.000000478233
Rust model of bus engine replacement (id=140628910063488) solved with poly in 81 iterations
At Rustโs original pure VFI would need tens of thousands of iterations, and the poly-algorithm solves the model in milliseconds โ which is what makes it possible to put the solver inside a likelihood on Thursday:
# original parameters from Rust 1987
m = zurcher()
polyset = {'sa_min':10,
'sa_max':100,
'switch_tol':0.0005,
}
ev,pk = m.solve_show(tol=1e-10,solver='poly',**polyset)
# time the solver
%timeit -n 5 -r 10 m.solve_poly(tol=1e-10,**polyset)
Rust model of bus engine replacement (id=140627722208128) solved with poly in 77 iterations
The slowest run took 420.18 times longer than the fastest. This could mean that an intermediate result is being cached.
1.5 s ยฑ 1.39 s per loop (mean ยฑ std. dev. of 10 runs, 5 loops each)
The switching rule compares the error ratio to . What happens if the tolerance
switch_tolis too loose, and what if it is too tight?NK is equivalent to policy iterations. Where in
solve_nkis the policy evaluation step, and where is the policy improvement?What in the model changes if the maintenance cost is made non-linear in mileage, and what in the code?
- Rust, J. (1987). Optimal Replacement of GMC Bus Engines: An Empirical Model of Harold Zurcher. Econometrica, 55(5), 999. 10.2307/1911259
- Keane, M. P., & Wolpin, K. I. (1997). The Career Decisions of Young Men. Journal of Political Economy, 105(3), 473โ522. 10.1086/262080
- Rust, J., & Phelan, C. (1997). How Social Security and Medicare Affect Retirement Behavior in a World of Incomplete Markets. Econometrica, 65(4), 781โ831. 10.2307/2171940
- Erdem, T., & Keane, M. P. (1996). Decision-Making under Uncertainty: Capturing Dynamic Brand Choice Processes in Turbulent Consumer Goods Markets. Marketing Science, 15(1), 1โ20. 10.1287/mksc.15.1.1
- Arcidiacono, P. (2004). Ability Sorting and the Returns to College Major. Journal of Econometrics, 121(1โ2), 343โ375. 10.1016/j.jeconom.2003.10.010
- Kennan, J., & Walker, J. R. (2011). The Effect of Expected Income on Individual Migration Decisions. Econometrica, 79(1), 211โ251. 10.3982/ECTA4657
- Eckstein, Z., & Wolpin, K. I. (1999). Why Youths Drop out of High School: The Impact of Preferences, Opportunities, and Abilities. Econometrica, 67(6), 1295โ1339. 10.1111/1468-0262.00081
- Hendel, I., & Nevo, A. (2006). Measuring the Implications of Sales and Consumer Inventory Behavior. Econometrica, 74(6), 1637โ1673. 10.1111/j.1468-0262.2006.00721.x
- Ackerberg, D. A. (2003). Advertising, Learning, and Consumer Choice in Experience Good Markets: An Empirical Examination. International Economic Review, 44(3), 1007โ1040. 10.1111/1468-2354.t01-2-00098
- Fosgerau, M., Frejinger, E., & Karlstrom, A. (2013). A Link Based Network Route Choice Model with Unrestricted Choice Set. Transportation Research Part B: Methodological, 56, 70โ80. 10.1016/j.trb.2013.07.012
- Francesconi, M. (2002). A Joint Dynamic Model of Fertility and Work of Married Women. Journal of Labor Economics, 20(2), 336โ380. 10.1086/338220
- Buchinsky, M., Gotlibovski, C., & Lifshitz, O. (2014). Residential Location, Work Location, and Labor Market Outcomes of Immigrants in Israel. Econometrica, 82(3), 995โ1054. 10.3982/ECTA10029
- Agarwal, N., Ashlagi, I., Rees, M. A., Somaini, P., & Waldinger, D. (2021). Equilibrium Allocations under Alternative Waitlist Designs: Evidence from Deceased Donor Kidneys. Econometrica, 89(1), 37โ76. 10.3982/ECTA17017
- Gillingham, K., Iskhakov, F., Munk-Nielsen, A., Rust, J., & Schjerning, B. (2022). Equilibrium Trade in Automobiles. Journal of Political Economy, 130(10), 2534โ2593. 10.1086/720463
- Iskhakov, F., Lee, J., Rust, J., Schjerning, B., & Seo, K. (2016). Constrained optimization approaches to estimation of structural models: Comment. Econometrica, 84(1), 365โ370. 10.3982/ECTA12605