Numerical Methods

The respy package contains several numerical components. We discuss each in turn.

Differentiation

Derivatives are approximated by forward finite differences and used by derivative-based optimization algorithms and the scaling procedure. Thus, the only possible setting for now is options_spec["derivatives"] = "forward-differences".

Integration

Integrals are approximated by Monte Carlo integration and occur in two different places:

  • The solution of the model requires the evaluation of \(E\max\). This integral is approximated using the number of random draws specified under options_spec["solution"]["draws"]. The same random draws are used for all integrals within the same period.

  • The estimation of the model requires the simulation of the choice probabilities to evaluate the sample likelihood. This integral is approximated using the number of random draws specified in the under options_spec["estimation"]["draws"]. The same random draws are used for all integrals within the same period.

Optimization

The estimation of the model involves the minimization of the simulated negative log-likelihood of the sample. The available optimizers depend on the version of the program. If you use the Python implementation, then the Powell (Powell, 1964) and BFGS (Norcedal and Wright, 2006) algorithms are available through their scipy implementations. For the Fortran implementation, we provide the BFGS and NEWUOA (Powell, 2004) algorithms. The algorithm to be used is specified under options_spec["optimizer"] and its default settings under options_spec["SCIPY-LBFGSB"] for example.

  • Preconditioning

    We implemented a diagonal scale-based preconditioner based on the gradient. To stabilize the routine, the user needs to specify a minimum value for the derivative approximation. The details are governed by the settings in options_spec["preconditioning"].

Function Approximation

We follow Keane and Wolpin (1994) and allow to alleviate the computational burden by calculating the \(E\max\) only at a subset of states each period and interpolating its value for the rest. We implement their proposed interpolation function:

\[\begin{align} E\max - \max E = \pi_0 + \sum^4_{j = 1} \pi_{1j} (\max E - \bar{V}_j) + \sum^4_{j = 1} \pi_{2j} \left(\max E - \bar{V}_j\right)^{\tfrac{1}{2}}. \end{align}\]

\(\bar{V}_j\) is shorthand for the expected value of the alternative-specific value function and \(\max E = \max_k\{\bar{V}_j\}\) is its maximum among the choices available to the agent. The \(\pi\)’s are time-varying as they are estimated by ordinary least squares each period. The subset of interpolation points for the interpolating function is chosen at random for each period. The number of interpolation points remains constant across all periods. The number of interpolation points is selected in the INTERPOLATION section of the initialization file.

Function Smoothing

We simulate the agents’ choice probabilities to evaluate the negative log-likelihood of the sample. With only a finite number of draws, there is always the risk of simulating a zero probability for an agent’s observed decision. So we implement the logit-smoothed accept-reject simulator as suggested by McFadden (1989). The scale parameter \(\tau\) is set in the options specification under options_spec["estimation"]["tau"].

Miscellaneous

We use the LAPACK library for all numerical algebra. The generation of pseudorandom numbers differs between the Python and Fortran implementations. While they are generated by the Mersenne Twister (Matsumoto and Nishimura, 1998) in Python, we rely on the George Marsaglia’s KISS generator (Marsaglia, 1968) in Fortran.