Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
26,391 characters · 11 sections · 14 citation commands
Jacobian-free Efficient Pseudo-Likelihood (EPL) Algorithm
Structural estimations of models with equilibrium constraints are prevalent in economics. Examples include dynamic discrete games with incomplete information, and dearing2024efficient recently proposed Efficient Pseudo Likelihood (EPL) estimator for estimating the model. The idea is to sequentially update structural parameters $\theta$ and nuisance parameters $Y$ (choice-specific value functions) by introducing Newton-like mappings. Unlike conditional choice probability estimators and Nested Pseudo Likelihood (NPL) estimators (aguirregabiria2007sequential), the EPL estimator is guaranteed to be efficient under regularity conditions. Moreover, MLE is a fixed point of the EPL iterations. Unlike the MPEC (Mathematical Programming with Equilibrium Constraints) approach (egesdal2015estimating), the algorithm does not require sophisticated software. Compared to the Nested-Fixed Point (NFXP) approach, the EPL is computationally fast and robust even when multiple equilibria exist.
One obstacle to implementing the EPL algorithm is the need to compute Jacobians of the equilibrium constraints concerning $Y$. First, analytically deriving the Jacobians is cumbersome and prone to coding mistakes, especially when considering complicated models.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {For instance, igami2017estimating introduced nonstationarity and sequential move assumptions in dynamic games.} Second, the Jacobian is an $|Y|\times|Y|$-dimensional matrix, the memory requirement and the computational cost can be high, especially when the dimension of the nuisance parameters is large.
The current study proposes a simple procedure to avoid the problem mentioned above. The idea is to combine numerical derivatives (for computing Jacobian-vector products) and the Krylov method (for solving linear equations).\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {Analogous idea has been used to apply Newton's method for solving nonlinear equations. See Section 3 of kelley2003solving.} With this strategy, we can completely avoid computing the Jacobians. I show good computational performance of the proposed method by numerical experiments.
Note that the proposed method is not necessarily as fast as the original EPL algorithm. Nevertheless, the proposed method does not require computing Jacobians, and need less storage than the original EPL. Moreover,the proposed method is still much faster than the NFXP. The proposed method would be worth trying even when practitioners want to experiment whether the EPL works well, before coding complicated Jacobians.
Efficient Pseudo-Likelihood (EPL) estimator is designed to estimate the following problem:
where $Q_{N}\left(\theta,Y\right)\equiv\frac{1}{N}\sum_{i=1}^{N}\ln Pr\left(w_{i}|\theta,Y\right)$. $i$ denotes a sample, and $w_{i}$ denotes an observed outcome. $Y$ is called a nuisance parameter. In the dynamic discrete game discussed in Section (ref), choice-specific value function $v$ corresponds to $Y$.
Here, we further assume linearity concerning parameters $\theta$: \[ G\left(\theta,Y\right)=H\left(Y\right)\theta+z\left(Y\right). \] Here, $H\left(Y\right)\equiv\left(h_{i}\left(Y\right)\right)_{i=1,\cdots,|\Theta|}$ is a $|\mathcal{Y}|\times|\Theta|$-dimensional matrix, and $z\left(Y\right)$ is a $|\mathcal{Y}|$-dimensional vector.
Algorithm (ref) shows the steps of the EPL algorithm.
In the EPL algorithm, Step 2(a) would be the largest obstacle to the implementation. In the step, we need to compute the values of $\left(\nabla_{Y}G\left(\widehat{\theta}_{k-1},\widehat{Y}_{k-1}\right)\right)^{-1}H\left(\widehat{Y}_{k-1}\right)$ and $\left(\nabla_{Y}G\left(\widehat{\theta}_{k-1},\widehat{Y}_{k-1}\right)\right)^{-1}z\left(\widehat{Y}_{k-1}\right)$. Coding $\nabla_{Y}G\left(\widehat{\theta}_{k-1},\widehat{Y}_{k-1}\right)$ is cumbersome, and prone to coding mistakes, especially for complicated models. In addition, computing $|Y|\times|Y|$-dimensional matrix $\nabla_{Y}G\left(\widehat{\theta}_{k-1},\widehat{Y}_{k-1}\right)$ requires much memory requirement and computation time when $|Y|$ is large.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {We might be able to mitigate the problem by utilizing sparsity. However, the strategy does not work for dense problems.}
The proposed procedure builds on the ideas of the Krylov method and numerical derivatives. We discuss them before describing the proposed method.
Suppose we want to solve $Ad=b$ for $d\in\mathbb{R}^{m}$, where $A\in\mathbb{R}^{m\times m}$, $b\in\mathbb{R}^{m}$. Algorithm (ref) shows the basic structure of the Krylov method.
As shown in Algorithm (ref), the matrix $A$ appears only in the form of Jacobian-vector product $A\bm{v}$, where $\bm{v}$ denotes a vector. It implies that we can run the Krylov method if a function (linear operator) $g:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ such that $g(\bm{v})=A\bm{v}$ is specified. GMRES algorithm, one kind of the Krylov method, is available in many programming languages (e.g., MATLAB, Scipy package in Python, IterativeSolvers package in Julia) with the feature.
Generally, by Taylor's theorem, twice differentiable function $f$ satisfies the following equations:
Hence, $\left\Vert \frac{f(x+\epsilon v)-f(x-\epsilon v)}{2\epsilon}-\left(\nabla_{x}f(x)\right)v\right\Vert \leq\epsilon^{2}O\left(\left\Vert p\right\Vert ^{3}\right)$ holds. By taking small $\epsilon$, we can approximate the value of $\left(\nabla_{x}f(x)\right)v$ by $\frac{f(x+\epsilon v)-f(x-\epsilon v)}{2\epsilon}$ well.
By combining the two ideas, we can compute $\left(\nabla_{Y}G\left(\widehat{\theta}_{k-1},\widehat{Y}_{k-1}\right)\right)^{-1}\left(H\left(\widehat{Y}_{k-1}\right)\right)$ and $\left(\nabla_{Y}G\left(\widehat{\theta}_{k-1},\widehat{Y}_{k-1}\right)\right)^{-1}\left(z\left(\widehat{Y}_{k-1}\right)\right)$ without explicitly computing $\nabla_{Y}G\left(\widehat{\theta}_{k-1},\widehat{Y}_{k-1}\right)$. Algorithm (ref) shows the proposed procedure.
Note that the convergence properties of the EPL algorithm do not change as long as $d^{H*}$ and $d^{z*}$ approximate $\left(\nabla_{Y}G\left(\widehat{\theta}_{k-1},\widehat{Y}_{k-1}\right)\right)^{-1}\left(H\left(\widehat{Y}_{k-1}\right)\right)$ and $\left(\nabla_{Y}G\left(\widehat{\theta}_{k-1},\widehat{Y}_{k-1}\right)\right)^{-1}\left(z\left(\widehat{Y}_{k-1}\right)\right)$ well. Concerning the numerical accuracy, see the discussion in Appendix (ref).
To show the performance of the proposed method, we consider a stationary dynamic discrete game of incomplete information as in aguirregabiria2007sequential, egesdal2015estimating, and dearing2024efficient.
Time is discrete, and indexed by $t=1,2,\cdots$. In each market, there are $|\mathcal{J}|$ firms, and they are indexed by $j\in\{1,2,\cdots,|\mathcal{J}|\}$. Given observable states $x_{t}$ and private information only observed to each firm $\epsilon_{t}^{j}$, each firm simultaneously chooses its action $a_{t}^{j}\in\mathcal{A}=\{0,1,\cdots,|\mathcal{A}|-1\}$.
Agents maximize the expected discounted utility $E\left\{ \sum_{s=0}^{\infty}\beta^{s-t}\left[\overline{u}^{j}\left(x_{s},a_{s}^{j},a_{s}^{-j};\theta_{u}\right)+\epsilon_{s}^{j}(a_{s}^{j})\right]\left|x_{t},\epsilon_{t}^{j}\right.\right\} $, where $\beta\in(0,1)$ denotes the discount factor. $\theta_{u}$ and $\theta_{f}$ represent parameters concerning utility and state transitions respectively, and let $\theta\equiv(\theta_{u},\theta_{f})$. Under the assumptions of conditional independence, independent private values, finite observed state space $\left(x_{t}\in\chi=\{1,2,\cdots,|\chi|\}\right)$, and Markov perfect equilibrium, choice-specific value function $v^{j}\in\mathbb{R}^{|\chi|\times|\mathcal{A}|}$ satisfies the following equation\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {See Lemma 1 of dearing2024efficient for the proof.}:
Here, we define a function $\Phi_{v}:\Theta\times\mathbb{R}^{|\mathcal{J}|\times|\chi|\times|\mathcal{A}|}\rightarrow\mathbb{R}^{|\mathcal{J}|\times|\chi|\times|\mathcal{A}|}$. $\Lambda^{j}(x,a^{j};v^{j})$ is the choice probability that agent $j$ chooses action $j$ in state $x$, conditional on having choice-specific value function $v^{j}$. $S(\cdot)$ is McFadden's surplus function. We also define $\Lambda^{-j}(v^{-j})\equiv\left(\Lambda^{1}(v^{1}),\cdots,\Lambda^{j-1}(v^{j-1}),\Lambda^{j+1}(v^{j+1}),\cdots,\Lambda^{|\mathcal{J}|}(v^{|\mathcal{J}|})\right)$, and:
Choice-specific value function $v$ is a solution of a nonlinear equation $G(\theta,v)\equiv v-\Phi_{v}(\theta,v)=0$ given $\theta$.
We consider the setting where firms make entry/exit decisions $(a^{j}=1,0)$. Concerning $\overline{u}^{j}$, we assume:
where $x_{t}\equiv\left(s_{t},\left\{ a_{t-1}^{j}\right\} _{j\in\mathcal{J}}\right)$ and $s_{t}\in\{1,\cdots,|\mathcal{S}|\}$ denotes the market size.
Then, $\overline{u}^{j}(x_{t},a_{t}^{j},a_{t}^{-j};\theta)=\exists h(x_{t},a_{t}^{j},a_{t}^{-j})^{\prime}\theta$ holds, and $G(\theta,v)$ can be represented as $G(\theta,v)=\exists H(v)\theta+\exists z(v)$ (linearity concerning $\theta$). We further assume $\epsilon$ follows i.i.d. type-I extreme value distribution.
The experiments are conducted on a laptop computer with the CPU AMD Ryzen 5 6600H 3.30 GHz, 16.0 GB of RAM, Windows 11 64-bit, and MATLAB 2022b. The experiments are based on the code building on the replication code of dearing2024efficient. As in dearing2024efficient, we choose parameter values $|\mathcal{J}|=5,|\mathcal{S}|=5,\theta_{FC,j}=-2+0.1\times j,\theta_{RS}=\theta_{RN}=\theta_{EC}=1.0$. Sample size $N$ is set to 1600. Concerning $\theta_{RN}$, we use $\theta_{RN}=4$. We assume the EPL iteration converges when $||\theta^{(n+1)}-\theta^{(n)}||_{\infty}$ and $||v^{(n+1)}-v^{(n)}||_{\infty}$ are less than (1E-2)$\slash K$, where $K$ denotes the number of parameters. Concerning the initial values, we initialize k-NPL with estimated semiparametric logit choice probabilities, and then initialize 1-EPL using the consistent parameter estimates and value function from the 1-NPL iteration.
As shown in the replication code of dearing2024efficient, we can use the following sparsity of the Jacobian $\nabla_{Y}G\left(\widehat{\theta}_{k-1},\widehat{Y}_{k-1}\right)$ when computing the matrix: $\frac{\partial\Phi^{j}(x_{t},a_{t}^{j};v^{j},v^{-j},\theta)}{\partial v^{l}\left(\widetilde{x}_{t},a_{t}^{l}\right)}=0\ \text{if}\ l\neq j\ \text{and}\ x\neq\widetilde{x}$. Though the Jacobian's total number of elements is $|\chi|^{2}|\mathcal{J}|^{2}|\mathcal{A}|^{2}$, we have to store only $|\chi|^{2}|\mathcal{J}||\mathcal{A}|^{2}+|\chi|\left(|\mathcal{J}|^{2}-|\mathcal{J}|\right)|\mathcal{A}|^{2}$ elements when utilizing the sparsity structure. In contrast, in the proposed Jacobian-free method, we only have to compute the $\left|\chi\right|\times\left|\chi\right|$-dimensional state transition matrix, which requires less memory.
To compute numerical derivatives in the proposed algorithm, I use $\epsilon=\frac{\sqrt[3]{\mathbf{u}}}{\max\left\{ ||v||_{\infty},\text{1E-8}\right\} }$, where $\mathbf{u}$ denotes the machine precision.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {Generally, $\left\Vert \frac{\widetilde{f}(x+\epsilon v)-\widetilde{f}(x-\epsilon v)}{2\epsilon}-\nabla f(x)^{T}v\right\Vert \leq\epsilon^{2}O\left(\left\Vert v\right\Vert ^{3}\right)+\frac{L_{f}}{\epsilon}$ holds for twice differentiable function $f$ when $\left\Vert \widetilde{f}(x)-f(x)\right\Vert \leq L_{f}\ \forall x$. The right hand side of the inequality is minimized when $\epsilon=\frac{L_{f}}{\sqrt[3]{O\left(\left\Vert v\right\Vert ^{3}\right)}}$. In the current setting, evaluating $G(\theta,v)$ incurs rounding errors, and $L_{f}$ would be around the machine precision $\mathbf{u}$. Hence, I choose the specification $\epsilon=\frac{\sqrt[3]{\mathbf{u}}}{\max\left\{ ||v||_{\infty},\text{1E-8}\right\} }$. Note that 1E-8 is introduced to avoid the division by 0 when $||v||_{\infty}$ equals 0.} As $\epsilon_{GMRES}$ (tolerance level in the GMRES iteration), I use 1E-5.
Table (ref) compares the proposed EPL algorithm with the original EPL algorithm. Here, $k$-EPL denotes the EPL estimator where the number of iterations is set to $k$. $\infty$-EPL denotes the estimator after the convergence of the iterations. “Anal.” denotes the EPL estimator using analytical Jacobians and mldivide function (non-iterative method to solve linear equations) in MATLAB. “Krylov” denotes the EPL estimator using analytical Jacobians and gmres function (iterative method to solve linear equations) in MATLAB. Finally, “JF” denotes the EPL estimator based on the proposed Jacobian-free approach. The results show that the difference between the EPL-JF and EPL-Analytical is mostly negligible. Concerning the computation time, EPL-JF is several times slower than EPL-Analytical and EPL-Krylov. Nevertheless, considering the benefit from avoiding coding Jacobians, EPL-JF would be attractive.
Note that EPL-JF requires much less computation time, compared to the NFXP algorithm using numerical derivatives in the outer-loop. The latter is also Jacobian-free, and we denote it as NFXP-JF. Table (ref) compares the performance,\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {Concerning the inner-loop of NFXP, I use the fixed-point iterations using the mapping $\Phi_{v}$. The tolerance level of the inner-loop iteration is set to 1E-12. To accelerate the convergence, I combine the Anderson acceleration method. Though I also experimented other fixed-point mappings and acceleration methods, the results were not largely different. Concerning the outer-loop optimization, the tolerance level is set to 1E-6, and I use fminunc function in MATLAB with central-difference numerical derivatives. In the current study, the inner-loop iterations are started from only one initial values of $v$. Generally, when we cannot rule out the possibility of multiple equilibria, we should start from several initial values. When we introduce such procedure, the computational cost of the NFXP algorithm would be much larger.} and we can see that the computational cost of EPL-JF is tens of times smaller than the NFXP.