The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
78,155 characters
Estimation and inference in models with multiple behavioural equilibria
\title{\vspace{-3em} Estimation and inference in models with multiple behavioural equilibria\thanks{We are grateful for the helpful suggestions and comments received at seminars at Erasmus University Rotterdam, Rijksuniversiteit Groningen, Collegio Sant’Anna (University of Pisa), and the CFE 2025 Conference (Birkbeck, University of London). We would like to thank Mariia Artemova, Giulio Bottazzi, John Cochrane, Pietro Dindo, Christian Francq, Rutger-Jan Lange, Sven Otto, Dario Palumbo, Andreas Pick and Dominik Wied for helpful comments and discussions. The second author acknowledges financial support under the National Recovery and
Resilience Plan (NRRP), Mission 4, Component 2, Investment 1.1, Call for tender No. 104 published
on 2.2.2022 by the Italian Ministry of University and Research (MUR), funded by the European Union
– NextGenerationEU– Project Title ”Market Learning and Robust Predictions (MALERP)” – CUP J53D23004190006. }}
\author{Alexander Mayer\footnote{{\it Corresponding author}, email: [email removed]}\\
{\it Erasmus University \& Tinbergen Institute}\\
Rotterdam, The Netherlands
\and
Davide Raggi\footnote{Email: [email removed]}\\
{\it Università Ca' Foscari}\\
Venezia, Italy}
\date{\today}
\maketitle
\thispagestyle{empty}
\vspace*{-.5cm}
\begin{abstract} We develop estimation and inference methods for a stylized macroeconomic model with potentially multiple behavioural equilibria, where agents form expectations using a constant-gain learning rule. We first show geometric ergodicity of the underlying process to study in a second step (strong) consistency and asymptotic normality of the nonlinear least squares estimator for the structural parameters. We propose inference procedures for the structural parameters and uniform confidence bands for the equilibria. When equilibrium solutions are repeated, mixed convergence rates and non-standard limit distributions emerge. Monte Carlo simulations and an empirical application illustrate the finite-sample performance of our methods.
\vspace{1em}\noindent{\bf Keywords:} adaptive learning, nonlinear least squares, ergodicity, consistency, asymptotic distribution, identification.
\end{abstract}
\section{Introduction}
The rational expectations hypothesis, proposed in the seminal work by \citet{Muth1961} and \citet{Lucas1972}, is a benchmark for the study of macroeconomic models, as it allows for a dynamic system in which agents' expectations are incorporated in a manner internally consistent with the model itself. Starting with \cite{marsar:89} and \citet{Sargent93,Sargent99}, however, the orthodox view of agents forming full-information rational expectations has been increasingly questioned; see, e.g., \citet{evans:01}, \citet{MankiwReis2002}, and \citet{Bernanke07}. With a particular focus on inflation expectations, a large body of recent empirical studies that exploit survey data corroborates the finding that agents deviate from the benchmark of full rationality; see, e.g., \citet{CoibionGorodnichenko2012}, \citet{BachmannBergSims15}, \cite{mn:2016}, \citet{CoibionGorodnichenkoRopele20}, \citet{coibion:20}, and \citet{BordaloGennaioliEtAl2020}.
Several deviations from the pure rational expectations paradigm have since been proposed. Notable examples include models with frictions, such as sticky information (e.g., \citealp{MankiwReis2002}) and rational inattention (e.g., \citealp{Sims2003}). Another popular route, drawing on bounded rationality and popularized by, among others, \cite{bray1982learning}, \cite{bray1986rational}, \cite{marsar:89}, and \cite{evans:01}, treats agents as statisticians who form expectations using simple, possibly misspecified, econometric models.
This is also the approach taken in this paper, which closely follows the specifications studied in \citet{HommesSorger1998}, \cite{Lansing2009}, \citet{hz14}, and \citet{HommesMavromatisOzdenZhu2023}. In particular, we consider a New Keynesian Phillips curve (NKPC) in which agents update their beliefs by recursively estimating the parameters of a first-order autoregression, AR(1), that they (incorrectly) perceive to be the correct model of the inflation process. As \citet{hz14} show theoretically, this framework provides a parsimonious account of several features of observed inflation dynamics.
Importantly, because agents remain even in the long run boundedly rational without correctly perceiving all stochastic aspects of the model, multiple so-called {\it behavioural equilibria} can coexist: self-confirming inflation paths that are fixed points of the perceived law of motion, rather than rational-expectations solutions of the full model. Recently, \citet{HommesMavromatisOzdenZhu2023} and \cite{davide:25} take variations of this model to the data and, using Bayesian techniques, show its usefulness both in-sample and out-of-sample.
As we discuss in detail below, however, because of the model’s nonlinear features and the presence of multiple behavioural equilibria, a priori little is known about the properties of (frequentist) estimators for the structural parameters in these models. A growing literature on the econometrics of models with adaptive learning shows that the statistical analysis is indeed a challenging task; see, e.g., \citet{chev:10}, \citet{adam:16}, \citet{chev:17}, \citet{chrismass:18,chrismass:19}, \citet{mayer:22,mayer:23}, \citet{christiano:24}, and \citet{mm:25}.
To close this gap, we {\it first} provide a thorough statistical analysis of a plain-vanilla NKPC model with constant-gain learning. Specifically, we $(i)$ establish geometric ergodicity of the dynamic system; $(ii)$ prove strong consistency and asymptotic normality of the nonlinear least-squares (NLS) estimator of the structural parameters and the associated behavioural equilibria; and $(iii)$ provide detailed guidance on conducting inference. Substantively, this delivers an estimable NKPC with learning and a characterization of the location of its behavioural equilibria. Methodologically, we develop a general strategy for establishing ergodicity and NLS inference in nonlinear learning models with potentially multiple equilibria, and highlight when and why standard limit theory breaks down.
In doing so, we draw on recent results for nonlinear time series models to derive the model’s probabilistic properties (\citealp{chot:19}) and to justify point estimation (\citealp{francq:24}). Moreover, several less-standard results emerge that connect to the broader econometric literature. For example, akin to \citet{hansen:1996}, we show how to conduct inference when certain parameters (viz., the learning gain) are not {\it jointly} identified. In this case, following \cite{saikkonen:95} and \cite{seo:11}, we can still verify $\sqrt n$-consistency of the remaining structural parameters. Our analysis of the (multiple) equilibria is also conceptually related to \citet{kasy:15}, who develops a nonparametric procedure for inference on the number of roots of an unknown function. In contrast, we work in a parametric setup where the equilibria are the roots of a polynomial function that is known up to a finite dimensional parameter. This allows us to study the asymptotic distribution of equilibrium locations, its number, and uniform inference on the complete function. When inference concerns the (multiple) equilibria, convergence can be slow, akin to the second-order identification mechanism in \citet{dovonon:18}. {\it Secondly}, we verify the theoretical findings in finite samples, using both simulated and real-world data. The empirical application to US data illustrates how our methods can be used to characterize uncertainty about equilibrium locations and learning dynamics in practice.
The remainder of this paper is organized as follows: Section \ref{sec:model} introduces the model. The time series properties of the resulting data generating process (dgp) are established in Section \ref{sec:prop}. Sections \ref{sec:cons} and \ref{sec:dis} discuss, respectively, consistency and asymptotic normality of the NLS estimator. The asymptotic distribution of the resulting estimator of the behavioural equilibria is derived in Section~\ref{sec:roots}, while Section~\ref{sec:inf} is devoted to various inference procedures. Finite sample evidence in the form of a Monte Carlo experiment and an estimation exercise using US data are presented in Sections~\ref{sec:mc} and \ref{sec:emp}, respectively, before Section~\ref{sec:conc} concludes. All proofs are delegated to the appendix.
\section{Model}\label{sec:model}
Consider a standard version of a forward-looking NKPC:
\begin{align}\label{eq:pc}
\pi_t = \delta \pi_{t+1}^e + \psi y_t + u_t,
\end{align}
where $t=1,2,\dots$ indexes time, $\pi_t$ denotes inflation, $\pi_{t+1}^e$ is the agents', not necessarily rational, expectation of future inflation, the driving variable $y_t$ is a proxy for real marginal costs (e.g. an output gap measure or real unit labour costs), and $u_t$ is an error term (see, e.g., \citealp{Woodford2003} for a theoretical discussion or \citealp{MavroeidisPlagborgMollerStock2014} for a more empirical treatment). Following, among others, \citet{Lansing2009} and \citet[Section 4]{hz14}, the driving variable is assumed to follow an AR(1) process,
\begin{align}\label{eq:AR1}
y_t = a + \rho y_{t-1} + \varepsilon_t,
\end{align}
with intercept $a$, persistence parameter $\rho$, and innovation $\varepsilon_t$. The parameters $\rho$ and $\delta$ are in the interval $(-1,1)$, whereas the variances of the shocks $\varepsilon_t$ and $u_t$ are given by $0<\sigma^2_\varepsilon < \infty$ and $0<\sigma^2_u < \infty$, respectively. The relevant economic parameter is $\psi \in {\mathbb R}$, usually called the slope, that measures how inflation responds to changes in economic `{slack}'.
Under the full information rational expectations (RE) hypothesis, agents have complete knowledge of the model structure (e.g. they know the functional form of Eq.~\eqref{eq:pc}) and its history and make optimal use of it by setting $\pi_{t+1}^e = \textnormal{\textsf{E}}[\pi_{t+1} \mid \mathcal{F}_t]$, where $\mathcal{F}_t \coloneqq \sigma(\{(u_i,\varepsilon_i): i \leq t\})$. It can easily be shown (see e.g. \citealp[Appendix B]{hz14}) that the RE equilibrium is given by a linear projection of $\pi_t$ on $(1,y_t)^{\textnormal{\textsf{T}}}$:
\begin{equation}\label{eq:ree}
\pi_t= \frac{a\delta \psi}{(1-\delta\rho)(1-\delta)}+ \frac{\psi}{(1-\delta\rho)}y_t + u_t.
\end{equation}
The RE paradigm, although providing an important benchmark, has been challenged over the years. One popular alternative, proposed by the macroeconomic learning literature, is to assume that agents themselves act like statisticians when forming their inflation expectations. In this literature, agents usually know the functional form of the RE equilibrium \eqref{eq:ree} but not its parameters, which they estimate; hence the term \emph{bounded rationality} (see, e.g., \citealp{evans:01} and \citealp{evan:20} for overviews).
Here, we deviate from complete rationality a bit further: Following \citet{hz14}, we assume that agents even lack knowledge of the functional form. Instead, expectations of tomorrow's inflation are formed according to a perceived law of motion (PLM) that incorrectly stipulates AR(1) dynamics for $\pi_t$:
\begin{align}\label{eq:PLM}
\pi_t = \alpha + \beta(\pi_{t-1}-\alpha) + e_t,
\end{align}
with $e_t$ denoting an error, while $(\alpha, \beta)^{\textnormal{\textsf{T}}}$ are parameters unknown by the agent, who has to estimate them. The PLM in Eq.~\eqref{eq:PLM} is thus an auxiliary model employed by the agent to update her beliefs about future inflation. Note that, similar to the discussion in \citet[Ch. 13]{evans:01}, the PLM is misspecified in the sense that it does not nest the RE equilibrium Eq.~\eqref{eq:ree}.
Two questions arise: (1) \emph{How do agents forecast?} and (2) \emph{What do agents actually forecast?} To address the first question, assume that expectations are formed at the end of period $t\!-\!1$ using a mean-squared-error-optimal two-period-ahead rule based only on information available at $t\!-\!1$\footnote{This mechanism--where forecasts are based on lagged rather than current information--is standard in the adaptive-learning literature (e.g., \citealp{DeGrauweMarkiewicz2013} or \citealp{LansingMa2017}) and is consistent with the timeline of events in \citet{elias2016}; see also \citet{adam2003learning} for a theoretical justification.}. This implies $\pi_{t+1}^e = \alpha + \beta^2 (\pi_{t-1}-\alpha)$ and, therefore, the implied actual law of motion (ALM) for inflation is
\begin{equation}\label{eq:iALM}
\pi_t = \delta[\alpha + \beta^2 (\pi_{t-1}-\alpha) ] + \psi y_t + u_t,
\end{equation}
which would obtain {\it if} agents knew the parameters $(\alpha,\beta)^{\textnormal{\textsf{T}}}$ governing their forecasting rule in Eq.~\eqref{eq:PLM}.
To answer the second question, \cite{hz14} introduce the notion of {\it behavioural learning equilibria}. With agents misspecifying the functional form (note that they disregard the regressor $y_t$ entirely), \cite{hz14} discipline behaviour from becoming overly irrational by imposing that ($a$) the unconditional mean and ($b$) the first-order autocorrelation of $\pi_t$ are correctly perceived. These two restrictions pin down the equilibrium values of $(\alpha,\beta)^{\textnormal{\textsf{T}}}$. In particular, \cite{hz14} equate the PLM and Eq.~\eqref{eq:iALM} on their unconditional moments, which yields that expected inflation under the ALM satisfies
\begin{equation}\label{eq:consE}
\alpha= \frac{\psi a}{(1-\rho)(1-\delta)},
\end{equation}
whereas the first-order autocorrelation $\beta$ solves the following quartic equation
\begin{equation}\label{eq:consAC1}
\beta = F(\beta;\lambda), \quad F(\beta;\lambda) \coloneqq \delta\beta^2+\frac{\psi^2\rho(1-\delta^2\beta^4)}{\psi^2(\delta\beta^2\rho+1)+(1-\rho^2)(1-\delta\beta^2\rho)(\sigma_u/\sigma_\varepsilon)^2},
\end{equation}
where $\lambda \coloneqq (\delta,\psi,\rho,\sigma_u^2,\sigma_\varepsilon^2)^{\textnormal{\textsf{T}}}$. Put differently, the unconditional mean and autocorrelation agents believe in must equal the unconditional mean and autocorrelation actually generated by the NKPC under those beliefs. It is worth noting that, in this case, expected inflation corresponds to the rational expectation case. Furthermore, Eq.~\eqref{eq:consAC1} might deliver multiple solutions, representing different expectational equilibria, depending on the value assigned to the structural NKPC parameters. For instance, the solution is unique when $\psi>0$ is large enough and $\sigma^2_u$ is small. As discussed by \cite{hz14}, multiple equilibria correspond to different long-run inflation persistence regimes, and policy or shocks can move the economy between them. Specifically, drawing on the findings presented in \citet{blanchard:16} and \citet{JorgensenLansing2025}, an equilibrium characterized by a low $\beta$ parameter implies strongly anchored expectations, a scenario that aligns with the original specification of the Phillips curve, namely, one defined in terms of inflation levels and the output gap. Conversely, a high $\beta$ parameter would be more consistent with an accelerationist formulation of the curve.
Agents typically do not know the parameters $(\alpha,\beta)^{\textnormal{\textsf{T}}}$ of their forecasting model Eq.~\eqref{eq:PLM}. Instead, and in line with the adaptive learning literature (\citealp{evans:01}), we assume agents employ a recursive updating scheme to estimate them; formally, a stochastic approximation algorithm (\citealp{ben:90}). In particular, agents are presumed to use a constant gain version of the sample autocorrelation learning (SAC) algorithm proposed by \cite{hz14} (see also \citealp{HommesMavromatisOzdenZhu2023}) to estimate the subjective parameters $(\alpha,\beta)^{\textnormal{\textsf{T}}}$:
\begin{equation}
\alpha_t = (1-\gamma)\alpha_{t-1}+\gamma \pi_t, \quad \beta_t = \frac{\sum_{i=1}^{t-1}(1-\gamma)^{t-1-i}(\pi_{i+1}-\alpha_i)(\pi_{i}-\alpha_i)}{\sum_{i=1}^t(1-\gamma)^{t-i}(\pi_i-\alpha_{i-1})^2}, \label{eq:learning}
\end{equation}
for some so-called `gain' or `learning' parameter $\gamma > 0$ and initial value $\pi_0=\alpha_0$. Note that $\gamma > 0$ implies that more weight is attached to recent data and that the recursive estimators do not converge in probability but, as discussed below, converge in distribution to non-degenerate weak limits. If, in contrast, $\gamma \rightarrow 0$, then the adaptive learning estimates can converge under mild conditions to one of the possible behavioural equilibria computed from Eq.~\eqref{eq:consAC1}, and in this respect the equilibrium is called locally stable. Constant-gain learning with $\gamma>0$ implies, in turn, persistent deviations from rationality and is often employed in empirical applications where it is the preferred choice due to its tracking ability in nonstationary environments (see, e.g. \citealp{milani:07} or \citealp{chev:10}).\footnote{\cite{markiewicz2014adaptive} find that constant-gain least squares typically yields lower mean squared forecast errors than decreasing-gain learning, particularly for macroeconomic series such as inflation, GDP growth, and unemployment.} Importantly, in their empirical estimation of an extended version of the NKPC discussed here and in \cite{hz14}, \cite{HommesMavromatisOzdenZhu2023} employ the constant gain specification of Eq.~\eqref{eq:learning}.
To sum up, under learning, the resulting data generating process (dgp) is thus given by Eqs. \eqref{eq:AR1} (for $y_t$) and \eqref{eq:learning} (for $\alpha_t$, $\beta_t$) in conjunction with the true ALM:
\begin{equation}\label{eq:ALM}
\pi_t = \delta[\alpha_{t-1} + \beta_{t-1}^2 (\pi_{t-1}-\alpha_{t-1}) ] + \psi y_t + u_t.
\end{equation}
The unknown parameters in this model are $\theta \coloneqq (\gamma,\delta,\psi)^{\textnormal{\textsf{T}}}$ together with $(a,\rho,\sigma_u,\sigma_\varepsilon)$, that determine the potentially multiple equilibria $\beta$, all of which we aim to estimate. This task is nontrivial because the dgp defined by \eqref{eq:AR1}, \eqref{eq:learning}, and \eqref{eq:ALM} is highly nonlinear. \footnote{Note that, although the dgp can be seen as an `{\it observation driven}' model in the sense of \cite{cox:81} (see also \citealp{blasques:24}), commonly used results in that literature to establish the time series properties do not apply in a straightforward way. For example, as discussed in more detail below, results based on stochastic recurrence equations like \citet[Thm 2.8]{strau:06} require a Lipschitz condition that fails here.} Nevertheless, as discussed in the next section, it can be cast as a non-linear state-space model, in the sense of \citet{meyn:12}, and shown to possess sufficient regularity needed for our later analysis.
\section{Probabilistic properties of the dgp}\label{sec:prop}
Before turning to estimation and inference for the structural parameters, we first establish the time series properties needed to ensure that this exercise is meaningful. In particular, we establish that the model, summarized by the ALM in Eq.~\eqref{eq:ALM}, viewed as a nonlinear Markov chain, is geometrically ergodic. This matters because, irrespective of initialization, this guarantees that the chain converges to a unique stationary distribution and insures the {\it invertibility} of the time series model (see, e.g. \citealp{strau:06} or \citealp{bla:18}). Moreover, this notion of ergodicity justifies the use of law of large numbers (LLN) and central limit theorem (CLT) results for sample averages and smooth functionals which are key ingredients for our later analysis of estimation and inference (for discussions of these concepts in the context of non-linear time series models see, e.g. \citealp{carrasco2002mixing} or \citealp{kristensen2005geometric}).
To set the stage, note that we can represent $\beta_t$ also recursively i.e.
\[
\beta_t = \beta_{t-1}+ \frac{\gamma}{r_t} ((\pi_t-\alpha_{t-1})(\pi_{t-1}-\alpha_{t-1})-\beta_{t-1}(\pi_t-\alpha_{t-1})^2),
\]
and
\[
r_t = r_{t-1}+\gamma ((\pi_t-\alpha_{t-1})^2-r_{t-1}).
\]
We follow a convention from \cite{hz14}, and assume that the agents' initial value equals the observed inflation; i.e. $\alpha_0 \stackrel{!}{=} \pi_0 \eqqcolon \operatorname*{\mathfrak{a}} \in \mathbb{R}$. Thus setting $r_1 = \gamma(\pi_1-\alpha_0)^2$, ensures $\beta_1 = 0$ for any $r_0, \beta_0 \in \mathbb{R}$. Now, using these two recursions, we define to that end the $5\times 1$ state vector
$s_t \coloneqq (\pi_t,\alpha_t,y_t,\beta_t,r_t)^{\textnormal{\textsf{T}}}$. In Appendix \ref{sec:A1} we show that the dynamic system $s_t$ can be viewed as a nonlinear Markov recursion $s_t = \mathcal{G}(s_{t-1},v_t),$ for a continuous mapping $\mathcal{G}(s,v)$ and innovation vector $v_t \coloneqq (u_t,\varepsilon_t)^{\textnormal{\textsf{T}}}$. It is worth noting that $s_t$ belongs to the class of nonlinear state space models, as defined in \citet[p. 33]{meyn:12} and that, by Assumption \ref{ass:density} below, $s_t$ is weak Feller (\citealp[Prop. 6.1.2]{meyn:12}).
We cannot, however, directly rely on off-the-shelf results typically used in the nonlinear time series literature to verify geometric ergodicity of $s_t$. For example, the otherwise very general results in \citet[Section~A.1]{tong1990non} do not apply because, among other reasons, the underlying skeleton difference equation fails to be Lipschitz, owing to the highly nonlinear dependence induced by $\beta_t$. For a similar reason, we cannot invoke \citet[Thm 2.8]{strau:06} or the ergodicity result from \citet[Lem 3]{ds:93} used, for instance, by \citet[Appendix D]{adam:16}; see Remark~\ref{rem:erg} for details. Instead, we derive the stochastic properties of the model from first principles, following the exposition in \citet{chot:19}.
In doing so, we impose the following assumption:
\begin{assumption}\label{ass:density}
The process $\varepsilon_t$ and $u_t$ are independently and identically distributed $(${\sf IID}$)$, with finite second moments, both admit lower semicontinous densities on ${\mathbb R}$, and $\varepsilon_t \perp u_t$.
\end{assumption}
While some mild deviations from the {\sf IID} assumption might in principle be possible, existence of lower semicontinous error-densities with full support is essential for both our ergodicity proof as well as parameter identification. Under Assumption~\ref{ass:density}, we can show that $s_t$ explores all meaningful regions of the state space ($\varphi$-irreducibility) and avoids periodic behaviour (aperiodicity). To obtain ergodicity we also need to preclude divergence to infinity. This can be ensured by a so-called {\it drift condition}, for which we require the following parameter constraints.
\begin{assumption}\label{ass:para}
$|\delta|<1$, $|\rho|<1$, $\psi \in \mathbb{R}$, and $\gamma\in(0, 1)$
\end{assumption}
The bounds $|\delta| <1$ and $|\rho|<1$ guarantee stability and ergodicity of the dynamic system, while $\psi$ can be left unrestricted for many of our core results. The gain is by construction positive, while $\gamma < 1$ is not restrictive as it covers the range of estimates of $\gamma$ typically found in empirical work (see, e.g. \citealp[Fig. 3]{berardi2017empirical}). We occasionally impose additional parameter restrictions to obtain further insights: for instance, $\rho > 0$ ensures existence of at least one equilibrium, $\delta \neq 0$ guarantees that $\gamma$ is {\it jointly} identified with $(\delta,\psi)$, while $\psi > 0$ is imposed to obtain valid uniform confidence bands. In line with economic theory, \citet[Section 4.2]{hz14} assume $\rho \in [0,1)$, $\psi>0$, and $\delta \in [0,1)$, where $\rho > 0$ is imposed to ensure existence of at least one equilibrium.
\begin{proposition}\label{prop:ergod}
If Assumptions~\textup{\ref{ass:density}} and \textup{\ref{ass:para}} hold, then the process $s_t$ is geometrically ergodic, with $\textnormal{\textsf{E}}[|r_t|]<\infty$ and $\textnormal{\textsf{E}}[\|z_t\|^2]<\infty$, where $z_t \coloneqq (\pi_t,\alpha_t,y_t)^{\textnormal{\textsf{T}}}$. Moreover, if, in addition, $\textnormal{\textsf{E}}[\|v_t\|^k] < \infty$, $v_t = (u_t,\varepsilon_t)^{\textnormal{\textsf{T}}}$, then $\textnormal{\textsf{E}}[\|z_t\|^k] < \infty$ and $\textnormal{\textsf{E}}[|r_t|^{k/2}]<\infty$ for $k>2$.
\end{proposition}
Intuitively, geometric ergodicity ensures exponentially fast forgetting of initial conditions, a crucial requirement for valid estimation and inference covered in the following section.
\begin{remark}\label{rem:abconv}
Another consequence of the preceding proposition is that the recursive estimators of the agent $($viz. $\alpha_t$ and $\beta_t)$ do not converge in probability for $\gamma > 0$ as $t\rightarrow \infty$; a known feature of constant-gain stochastic approximation algorithms $($see, e.g. \textup{\citealp{ben:90})}. Instead, $(\alpha_t,\beta_t)$ converge in distribution to a stationary law. The limiting law as $t \rightarrow \infty$ of $\alpha_t \rightarrow_d \alpha^\star$ is centred at the RE solution with finite variance, i.e.
\[
{\textnormal{\textsf{E}}}[\alpha^\star] = \frac{a\psi}{(1-\rho)(1-\delta)}, \quad \textnormal{\textsf{var}}[\alpha^\star] = \frac{\gamma}{2-\gamma}\sum_{i=-\infty}^\infty \textnormal{\textsf{cov}}[\pi,\pi_{-i}](1-\gamma)^{|i|}.
\]
A similar closed-form characterization of the first two moments of the stationary law of the recursive sample autocorrelation $\beta_t$ is not feasible due its nonlinear form in Eq.~\eqref{eq:learning} defined as a ratio.\footnote{However, shifting to the first auto\textup{covariance} $\omega_t \coloneqq r_t\beta_t$, yields tractable moments of the stationary law $\omega_t \rightarrow_d \omega^\star$: i.e., we get $\textnormal{\textsf{E}}[\omega^\star] = \textnormal{\textsf{cov}}[x,x_{-1}]$ and $\textnormal{\textsf{var}}[\omega^\star] = \frac{\gamma}{2-\gamma}\sum_{i=-\infty}^\infty \textnormal{\textsf{cov}}[w,w_{-i}](1-\gamma)^{|i|}$, where $w_t \coloneqq x_tx_{t-1}$.} Following \textup{\citet[Ch. 4]{ben:90}}, more could be said using a so-called `small-$\gamma$' Gaussian approximation: i.e., both $\gamma^{-1/2}(\alpha_t-\textnormal{\textsf{E}}[\alpha^\star])$ and $\gamma^{-1/2}(\beta_t-\textnormal{\textsf{corr}}[x,x_{-1}])$, $x_t \coloneqq \pi_t-\alpha_{t-1}$, converge in distribution to mean-zero Gaussian variates as
$\gamma \rightarrow 0_+$ such that $\gamma t \rightarrow \infty$; see also \textup{\citet[Section 7.5]{evans:01}}. We will, however, maintain throughout the assumption of a constant gain viewing $\gamma>0$ as a parameter to be estimated. For robust procedures covering the case $\gamma \rightarrow 0$, see \textup{\cite{chev:10}}.
\end{remark}
\section{Estimation and inference}\label{sec:est}
Observing a sample $\{\pi_t,y_t\}_{t=1}^n$ of length $n$, we estimate the $3 \times 1$ parameter vector $\theta \coloneqq (\gamma,\delta,\psi)^{\textnormal{\textsf{T}}}$, using the nonlinear least squares (NLS) estimator $\theta_n \coloneqq (\theta_{\gamma,n},\theta_{\delta,n},\theta_{\psi,n})^{\textnormal{\textsf{T}}}$, which minimizes, over a suitable parameter space $\Theta \subset {\mathbb R}^3$, the sample objective
\begin{align}
Q_n(\theta,\operatorname*{\mathfrak{a}}) \coloneqq \sum_{t=1}^n(\pi_{t} - f_t(\theta,\operatorname*{\mathfrak{a}}))^2,
\end{align}
where, for a given initial value $\operatorname*{\mathfrak{a}}$, $f_t(\cdot,\operatorname*{\mathfrak{a}})$ is the nonlinear regression function
\begin{align} \label{eq:regfun}
f_t(\theta,\operatorname*{\mathfrak{a}}) \coloneqq \delta[\alpha_{t-1}(\gamma,\operatorname*{\mathfrak{a}})+\beta_{t-1}^2(\gamma,\operatorname*{\mathfrak{a}})(\pi_{t-1}-\alpha_{t-1}(\gamma,\operatorname*{\mathfrak{a}}))]+\psi y_t.
\end{align}
Importantly, by Proposition~\ref{prop:ergod} the choice of initial condition is asymptotically irrelevant: in terms of \cite{bougerol:93}, \cite{strau:06}, or \cite{blasques:22}, the {\it filter} forgets $\operatorname*{\mathfrak{a}}$ at a geometric rate and we suppress in the following the dependence on $\operatorname*{\mathfrak{a}}$.
Given $\theta_n$, we can estimate $\sigma_u^2$ using $\hat\sigma_u^2 \coloneqq \frac1{n}Q_n(\theta_n)$. Moreover, let $\rho_n$ and $\hat\sigma_e^2$ denote the OLS estimators of the AR(1) process in Eq.~\eqref{eq:AR1} and its innovation variance $\sigma_e^2$, respectively. We can then define for $$\lambda_n \coloneqq (\theta_{\delta,n},\theta_{\psi,n},\rho_n,\hat\sigma_u^2,\hat\sigma_\varepsilon^2)^{\textnormal{\textsf{T}}},$$ the estimator(s) of the (potentially multiple) equilibria
\begin{align}\label{eq:betahat}
\vartheta_n = F(\vartheta_n;\lambda_n),
\end{align}
implicitly given as the solution of Eq.~\eqref{eq:consAC1}.
\begin{remark}\label{rem:prof} As in \textup{\cite{hansen:17} (}see also \textup{\citealp{mm:25})}, we can reduce the numerical complexity by profiling with respect to $\gamma$. That is, set
$(\delta_n,\psi_n)^{\textnormal{\textsf{T}}} \coloneqq (\delta_n(\gamma_n),\psi_n(\gamma_n))^{\textnormal{\textsf{T}}}$ and let $\gamma_n$ the estimator that minimises the profiled objective
$$Q^\star_n(\gamma) \coloneqq \sum_{t = 1}^{n}(\pi_{t} - \delta_n(\gamma) h_{t-1}(\gamma)-\psi_n(\gamma)y_{t})^2,$$
where, for a given $\gamma$, the $2\times 1$ OLS estimator is given by
$$(\delta_n(\gamma),\psi_n(\gamma))^{\textnormal{\textsf{T}}} \coloneqq \left[\sum_{t=1}^n (h_{t-1}(\gamma),y_t)^{\textnormal{\textsf{T}}}(h_{t-1}(\gamma),y_t)\right]^{-1}\sum_{t=1}^n(h_{t-1}(\gamma),y_t)^{\textnormal{\textsf{T}}}\pi_t.$$
\end{remark}
\subsection{Consistency of the NLS estimator}\label{sec:cons}
Consistency of the NLS estimator requires identification of the true parameter vector $\theta_0$. Typically (e.g.\ \citealp[Assumption~($b$)]{jennrich:1969}) one shows that the population objective
\(Q(\theta)\coloneqq\textnormal{\textsf{E}}[Q_n(\theta)/n]\) is uniquely minimized at \(\theta_0\).
In our setting this route is impractical: the nonlinearity induced by \(\beta_t\)—itself a ratio of
discounted partial sums that depend on \(\theta\)—renders \(Q(\theta)\) analytically intractable.
Instead, we adapt the strong-consistency argument in \citet[pp.~1445--1447]{francq:24}, which
extends \citet[Theorem~1]{francq:04} and does not require an explicit evaluation of the population
objective.
To illustrate, note that
$$\pi_t = a\psi+\delta h_{t-1}(\gamma)+\rho\psi y_{t-1}+w_t, \quad h_t(\gamma) \coloneqq \alpha_t(\gamma)+\beta_t(\gamma)^2(\pi_t-\alpha_t(\gamma)),$$
where \(w_t \coloneqq \psi\varepsilon_t+u_t\) drives the conditional law of \(s_t\mid\mathcal F_{t-1}\). A key ingredient for identification is to show that $h_t(\gamma)$ is well separated from $h_t(\gamma_0)$ whenever $\gamma \neq \gamma_0$. While the nonlinearity in $\beta_t(\gamma)$ complicates this task, on tail events \(\{|w_t|\ge c\}\) (which have strictly positive probability by the full-support
part of Assumption~\ref{ass:density}), the easier-to-handle $\alpha_t(\gamma)$ dominates. This yields the following separation result:
\begin{lemma}\label{lem:tail-sep}
If Assumption~\textup{\ref{ass:density}} and \textup{\ref{ass:para}} holds, then for any \(\gamma\neq\gamma_0\) there exists
a constant \(c>0\) such that
$\mathbb{P}\{|h_t(\gamma)-h_t(\gamma_0)|>c\} > 0.$
\end{lemma}
This result, combined with compactness of $\Theta$, yields strong consistency via the arguments of
\citet{francq:24}. We therefore impose the additional restriction of the parameter space.
\begin{assumption}\label{ass:para1}
$|\rho|<1$ and $\theta \in \Theta \coloneqq \Gamma \times \Delta \times \Psi,$ where $$\Gamma \coloneqq [\ubar\gamma,\bar\gamma], \quad \Delta \coloneqq [-\bar\delta,-\ubar\delta] \cup [\ubar\delta,\bar\delta], \quad \Psi\coloneqq [\ubar\psi,\bar\psi]$$for $0<\ubar\gamma<\bar\gamma<1$, $0<\ubar\delta<\bar\delta<1$, $-\infty < \ubar\psi < \bar\psi < \infty$.
\end{assumption}
We are now ready to state the first core result:
\begin{proposition}\label{prop:const} Suppose Assumptions~\textup{\ref{ass:density}} and \textup{\ref{ass:para1}} hold, then $\theta_n \rightarrow_{a.s.} \theta_0$ as $n\rightarrow \infty.$
\end{proposition}
It is immediate that $\delta = 0$ implies that $\gamma$ is no longer {\it jointly} identified with $(\delta,\psi)$. However, we can still show that the remaining components of the NLS estimator remain (weakly) $\sqrt n$-consistent. In particular, following \cite{saikkonen:95} and \cite{seo:11}, the idea is to show that under $\delta_0=0$
\begin{align*}
D_n(\theta)\coloneqq Q_n(\theta)-Q_n(\theta_0)
= \,& \eta^{\textnormal{\textsf{T}}} \sum_{t=1}^n z_t(\gamma)z_t(\gamma)^{\textnormal{\textsf{T}}} \eta
- 2\,\eta^{\textnormal{\textsf{T}}} \sum_{t=1}^n z_t(\gamma)u_t > 0,\\
\eta \coloneqq \,&(\delta,\psi-\psi_0)^{\textnormal{\textsf{T}}}, \; z_t(\gamma)\coloneqq (h_{t-1}(\gamma),y_t)^{\textnormal{\textsf{T}}},
\end{align*}
uniformly outside any shrinking neighbourhood of $\eta_0$ whose radius goes to zero slower than $1/\sqrt{n}$. This, in turn, forces the NLS minimizer to lie within an $1/\sqrt{n}$ neighbourhood of $\eta_0$. Our proof builds, amongst others, on a uniform LLN and CLT for $Z_n(\gamma) \coloneqq \frac1{n}\sum_{t=1}^n z_t(\gamma)z_t(\gamma)^{\textnormal{\textsf{T}}}$ and $S_n(\gamma) \coloneqq \frac1{\sqrt n}\sum_{t=1}^n z_t(\gamma) u_t$, respectively, which--following \citet[Lem 1]{andrews:92} and \citet[Thm 2]{hansen:1996b}--can be shown by establishing stochastic Lipschitz bounds. Technically, these Lipschitz conditions are a result of the following uniform moment bounds that will be used throughout:
\begin{lemma}\label{lem:mom}
Let $k \geq 1$ and denote by $\alpha_t^{(m)}(\gamma) \coloneqq \frac{{\sf d}^m}{{\sf d} \gamma^m}\alpha_t(\gamma)$ $($similarly for $r_t$ and $\beta_t$$)$ the $m$-th derivative with respect to $\gamma$.
\begin{enumerate}
\item[$(i)$] If $\textnormal{\textsf{E}}[\|v_t\|^k] < \infty$, then
$\textnormal{\textsf{E}}[\operatorname*{\textnormal{\textsf{sup}}}\limits_{\gamma \in \Gamma}|\alpha_t^{(m)}(\gamma)|^k] < \infty$, $m \in \{0,1,2,3\}.$
\item[$(ii)$] If $\textnormal{\textsf{E}}[\|v_t\|^{2k}] < \infty$, then
$\textnormal{\textsf{E}}[\operatorname*{\textnormal{\textsf{sup}}}\limits_{\gamma \in \Gamma}|r_t^{(m)}(\gamma)|^k] < \infty,$ $m \in \{0,1,2,3\}.$
\item[$(iii)$] If $\textnormal{\textsf{E}}[\|v_t\|^{4k}] < \infty$ and $\operatorname*{\textnormal{\textsf{inf}}}\limits_{\gamma \in \Gamma}r_t(\gamma) \geq r > 0$ $a.s.$, then
$\textnormal{\textsf{E}}[\operatorname*{\textnormal{\textsf{sup}}}\limits_{\gamma \in \Gamma}|\beta_t^{(m)}(\gamma)|^k] < \infty$, $m \in \{1,2\},$
and, if in addition, $\textnormal{\textsf{E}}[\|v_t\|^{6k}] < \infty$, then $\textnormal{\textsf{E}}[\operatorname*{\textnormal{\textsf{sup}}}\limits_{\gamma \in \Gamma}|\beta_t^{(3)}(\gamma)|^k] < \infty.$
\end{enumerate}
\end{lemma}
Intuitively, bounding moments of the derivatives of $r_t$ and $\beta_t$ requires more stringent moment conditions on the errors $v_t = (u_t,\varepsilon_t)^{\textnormal{\textsf{T}}}$ because $r_t$ and $\beta_t$ are functions of sample second moments and thus already involve squared data.
This leads us to impose the following restriction on the errors:
\begin{assumption}\label{ass:eight} $\textnormal{\textsf{E}}[\|v_t\|^{8}] < \infty$.
\end{assumption}
\begin{corollary}\label{cor:rate}
Suppose Assumptions~\textup{\ref{ass:density}}, \textup{\ref{ass:para}}, and \textup{\ref{ass:eight}} hold. If, in addition $\delta_0 = 0$ and $\gamma \in \Gamma$, then $\sqrt{n}(\theta_{n,\delta},(\theta_{n,\psi}-\psi_0))^{\textnormal{\textsf{T}}} = O_p(1).$
\end{corollary}
As will be shown in the following section, if $\delta_0 \neq 0$, then the joint NLS estimator (weakly) $\sqrt n$-consistent for the full parameter vector $(\gamma,\delta,\psi)^{\textnormal{\textsf{T}}}$ under a slightly strengthened set of assumptions.
\subsection{Asymptotic normality of the NLS estimator}\label{sec:dis}
Having established the consistency of $\theta_n$, we turn to the question of its limiting distribution. To aid the discussion, define the gradient vector and the Hessian matrix of the nonlinear regression function $f_t(\theta) = \delta h_{t-1}(\gamma)+\psi y_{t}$ with respect to $\theta = (\delta,\gamma,\psi)^{\textnormal{\textsf{T}}}$ by
\[
\dot f_t(\theta) = (h_{t-1}(\gamma),\delta_0 \dot{h}_{t-1}(\gamma),y_{t})^{\textnormal{\textsf{T}}} \quad \text{and} \quad H_{t}(\theta) \coloneqq \begin{bmatrix}
0 & \dot{h}_{t-1}(\gamma) & 0 \\
\dot{h}_{t-1}(\gamma) & \delta \ddot{h}_{t-1}(\gamma) & 0 \\
0 & 0 & 0
\end{bmatrix},
\]
where
$
\dot{h}_t(\gamma) \coloneqq (1-\beta_t(\gamma)^2)\dot{\alpha}_t(\gamma)+2\beta_t(\gamma)\dot{\beta}_t(\gamma)(\pi_t-\alpha_t(\gamma))
$
and
\begin{align*}
\ddot{h}_t(\gamma) \coloneqq & (1-\beta_t(\gamma)^2)\ddot{\alpha}_t(\gamma)\\
&\quad+2(\dot{\beta}_t(\gamma)^2+\ \beta_t(\gamma)\ddot{\beta}_t(\gamma))(\pi_t(\gamma)-\alpha_t(\gamma))- 4\beta_t(\gamma)\dot{\beta}_t(\gamma)\dot{\alpha}_t(\gamma).
\end{align*}
The notation $\dot \alpha$ and $\ddot \alpha$ ($\dot \beta$ and $\ddot \beta$) is used to indicate the first and second derivative wrt $\gamma$ of $\alpha$ ($\beta$). Using the definition of $\theta_n$ in conjunction with the mean-value theorem, ensures that there exists a $\bar\theta_n$ on the line segmenting connecting $\theta_n$ and $\theta_0$ such that
\[
\sqrt{n}(\theta_n-\theta_0) = - M_n(\bar\theta_n)^{-1}\frac1{\sqrt n}\sum_{t=1}^n \dot f_tu_t, \quad \dot f_t \coloneqq \dot f_t(\theta_0),
\]
where
\(
M_n(\theta) \coloneqq -\frac1{n}\sum_{t=1}^n \dot f_t(\theta)\dot f_t(\theta)^{\textnormal{\textsf{T}}} + \frac1{n}\sum_{t=1}^n(\pi_t-f_t(\theta))H_t(\theta).
\)
Given the consistency of $\theta_n$, standard arguments (see, e.g. \citealp{amemiya1985advanced}) yield limiting normality if the following holds true:
\begin{itemize}
\item[($a$)] $A \coloneqq \textnormal{\textsf{E}}[\dot f_t \dot f_t^{\textnormal{\textsf{T}}}]$ is a finite positive definite $3\times 3$ matrix;
\item[($b$)] $M_n(\bar\theta_n) \rightarrow_p -A$ for any $\bar\theta_n$ such that $|\bar\theta_n-\theta_0| \leq \varepsilon_n \rightarrow_{a.s.} 0$;
\item[($c$)] $\frac1{\sqrt n}\sum_{t=1}^n \dot f_tu_t \rightarrow_d {\cal N}(0,A\sigma_u^2)$.
\end{itemize}
Often, ($a$) can be ensured by direct evaluation of the expectation $A$. Here, due to the highly nonlinear nature of the $\dot f_t(\theta)$, deriving $A$ in closed form is not feasible. Still, we are able to verify ($a$) by a delicate tail argument, which, as an immediate consequence of Lemma~\ref{lem:pd} and the CLT for homoskedastic martingale difference sequences, yields part ($c$). Moreover, as argued in \citet[Lemma 1]{andrews:92}, part ($b$) follows from a stochastic Lipschitz bound. These results are summarized by the following lemma:
\begin{lemma}\label{lem:pd} Let $v_t = (u_t, \varepsilon_t)^{\textnormal{\textsf{T}}}$.
\textup{($a$)} If $\textnormal{\textsf{E}}[\|v_t\|^4] < \infty$, then $A$ is a finite and positive definite $3 \times 3$ matrix.
\textup{($b$)} If $\textnormal{\textsf{E}}[\|v_t\|^{8}]<\infty$, then $\|M_n(\theta_2)-M_n(\theta_1)\| \leq B_n \|\theta_1-\theta_2\|$, with $B_n = O_p(1)$ for any $\theta_1,\theta_2 \in \Theta$.
\end{lemma}
In view of Lemma~\ref{lem:mom}, it is interesting to note that, although part ($b$) involves bounds on expected moments of up to the third derivatives of the agents' recursive estimates, we do not require more than the eight error-moments already imposed by Assumption~\ref{ass:eight}.
We are now ready to derive the limiting distribution of the estimator:
\begin{proposition}\label{prop:norm} If Assumptions~\textup{\ref{ass:density}}, \textup{\ref{ass:para1}}, and \textup{\ref{ass:eight}} hold and $\theta_0 =(\gamma_0,\delta_0,\psi_0)^{\textnormal{\textsf{T}}} \in {\sf int}(\Theta)$, then $${\sqrt n}(\theta_n-\theta_0) \rightarrow_d {\mathcal N}(0, \Sigma), \quad \Sigma \coloneqq \sigma_u ^2 A^{-1} \succ 0.$$
\end{proposition}
As is common for extremum estimators (see, e.g. \citealp{newmc:94}), we restrict the true value $\theta_0$ to fall within the interior of $\Theta$, excluding boundary cases such like $\delta = 0$, which--as already discussed--require additional attention.
\subsection{Asymptotic distribution of the roots and the number of roots}\label{sec:roots}
Recall from Eq.~\eqref{eq:betahat} that an equilibrium is defined via $\vartheta_n = F(\vartheta_n;\lambda_n)$ for $\lambda_n = (\theta_{\delta,n},\theta_{\psi,n},\rho_n,\hat\sigma_u ^2,\hat\sigma_\varepsilon^2)^{\textnormal{\textsf{T}}}$, where $\theta_n =(\theta_{\gamma,n},\theta_{\delta,n}, \theta_{\psi,n})^{\textnormal{\textsf{T}}}$, while $\rho_n$ and $\hat\sigma_\varepsilon^2$ denote the OLS estimator of $\rho$ and the associated error variance $\sigma_\varepsilon^2$, respectively. In order to study the properties of our estimator of the equilibria, we thus need the following intermediate result:
\begin{corollary} Suppose the conditions of Proposition~\textup{\ref{prop:norm}} are satisfied, then $$\sqrt{n}(\lambda_n - \lambda_0) \rightarrow_d {\mathcal N}(0,\Omega), \quad \Omega \succ 0.$$
\end{corollary}
If $u_t$ is symmetric, $\Omega$ has a simple block-diagonal structure:
$$
\Omega = \begin{bmatrix}
\Sigma_{23} & 0 \\ 0 & D
\end{bmatrix},
$$
where $\Sigma_{23}$ is the $2 \times 2$ $(2,3)$-block of $\Sigma$
and $D \coloneqq \text{diag}(1-\rho^2,\textnormal{\textsf{E}}[u^4]-\sigma_u^4,\textnormal{\textsf{E}}[\varepsilon^4]-\sigma_\varepsilon^4)$.
Let $r(\lambda)$ denote the number of roots of $G(\beta;\lambda)=\beta - F(\beta;\lambda)$ for a given $\lambda$, and define $$r_0 \coloneqq r(\lambda_0),\quad r_n \coloneqq r(\lambda_n),\quad G(\beta) \coloneqq G(\beta;\lambda_0),\quad G_n(\beta) \coloneqq G(\beta;\lambda_n).$$ By Assumption~\ref{ass:para}, $G(\beta)$ is continuous on $[0,1]$, and, if in addition $\rho>0$ and $\psi \neq 0$, satisfies $G(0)<0<G(1)$. This implies $r_0 \in \{1,2,3\}$ and that the number of simple roots is odd. In particular, if $r_0=2$, then there is one simple root and one repeated root, at which $G(\cdot)$ does not change sign; this case is illustrated by the solid black line in Figure~\ref{fig:repeated}. Moreover, when $r_0=2$, the event $\{r_n=r_0\}$ has probability zero; in terms of Figure~\ref{fig:repeated}, one observes with equal probability either the dashed curve ($r_n=1$) or the dotted curve ($r_n=3$) for the realised $G_n(\beta)$. Because the derivative of $G(\beta)$ at the repeated root is zero, the equilibria are not identified at first order but only at second order, echoing the discussion in \cite{dovonon:18}. This leads to a slower convergence rate and is summarised by the following corollary.
\begin{figure}[tb]
\centering
\begin{minipage}[c]{.7\textwidth}
\begin{tikzpicture}
\begin{axis}[
width=\linewidth,
height=6cm,
xlabel={$\beta$},
ylabel={$G(\beta)$},
legend pos=south east,
legend cell align=right,
legend style={draw=none, fill=none,font=\footnotesize},
table/col sep=space,
grid style={line width=.1pt, draw=gray!20},
major grid style={line width=.2pt, draw=gray!40},
scaled y ticks=false,
yticklabel style={/pgf/number format/fixed, /pgf/number format/precision=4},
enlarge x limits=false,
xticklabels={},ytick={0},
yticklabels={0},
xlabel near ticks,
ylabel near ticks,
label style={font=\footnotesize},
tick label style={font=\scriptsize},
xlabel style={yshift=2pt},
ylabel style={xshift=2pt}
]
\draw[red, line width=.5pt]
(axis cs:\pgfkeysvalueof{/pgfplots/xmin}, 0)
-- (axis cs:\pgfkeysvalueof{/pgfplots/xmax}, 0);
\draw[gray, line width=.5pt]
(axis cs:0.7599, \pgfkeysvalueof{/pgfplots/ymin})
-- (axis cs:0.7599, \pgfkeysvalueof{/pgfplots/ymax});
\draw[gray, line width=.5pt]
(axis cs:0.6726, \pgfkeysvalueof{/pgfplots/ymin})
-- (axis cs:0.6726, \pgfkeysvalueof{/pgfplots/ymax});
\addplot+[smooth, dashed, no marks, thick, black] table[x=ggrid, y=g1] {dat.txt}; \addlegendentry{$r_n=1$}
\addplot+[smooth, no marks, ultra thick, black] table[x=ggrid, y=g] {dat.txt}; \addlegendentry{$r_0=2$}
\addplot+[smooth, dotted, no marks, thick, black] table[x=ggrid, y=g2] {dat.txt}; \addlegendentry{$r_n=3$}
\end{axis}
\end{tikzpicture}
\end{minipage}
\begin{minipage}[c]{.275\textwidth}\vspace*{-.35cm}
\caption{The solid black line is the true $G(\beta)$ at $r_0 = 2$ with two roots indicated by the two grey vertical lines. The dashed and the dotted lines indicate a realized $G_n(\beta)$ conditioning on $r_n = 1$ and $r_n = 3$, respectively.}\label{fig:repeated}
\end{minipage}
\end{figure}
To set the stage, recall that $\sqrt{n}(\lambda_n-\lambda_0) \rightarrow_d \mathbb{Z}\sim {\mathcal N}(0,\Omega)$, let $\vartheta_{i,n}$, $1 \leq i \leq 3$, denote the potential roots of $\beta \mapsto G_n(\beta)$ and impose the following:
\begin{assumption}\label{ass:para2} Assumption~\textup{\ref{ass:para1}} holds with $\rho \in (0,1)$ and $\Theta \coloneqq \Gamma \times \Delta \times \Psi_{\neq 0}$, where $$\Psi_{\neq 0} \coloneqq [-\bar\psi,-\ubar\psi] \cup [\ubar\psi,\bar\psi],\quad 0<\ubar\psi<\bar\psi<\infty.$$
\end{assumption}
Assumption~\ref{ass:para2} introduces the additional restriction $\psi \neq 0$. Since $\|F_\lambda(\beta;\lambda)\| = \beta^2$ if the true $\psi_0 = 0$, this is important to ensure that $F_\lambda(\beta;\lambda)$ is uniformly bounded away from zero at $\lambda = \lambda_0$; a crucial condition for the later analysis. The sign restriction $\rho > 0$ together with $\psi \neq 0$ ensures $G(0) < 0 < G(1)$ (e.g. existence of at least one root).
\begin{corollary}\label{cor:roots} Suppose Assumptions~\textup{\ref{ass:density}}, \textup{\ref{ass:para2}}, and \textup{\ref{ass:eight}} hold.
\begin{itemize}
\item[\textnormal{1)}]
If $r_0 \in \{1,3\}$, $r_n \rightarrow_p r_0$ and
\[
\sqrt{n}(\vartheta_{i,n}-\beta_{i,0}) \rightarrow_d J_i^{\textnormal{\textsf{T}}}{\mathbb Z},
\quad
J_i\coloneqq \frac{F_\lambda(\beta_{i,0};\lambda_0)}{\,1-F_\beta(\beta_{i,0};\lambda_0)\,}, \quad i =1,\dots, r_0.
\]
\item[\textnormal{2)}] If $r_0 = 2$, $r_n \rightarrow_d {\sf Bern}(\frac1{2})$ on $\{1,3\}$, i.e. $r_n \not\to_p r_0$.
\begin{itemize}
\item[\textnormal{($a$)}] If $r_n = 1$, $\sqrt{n}(\vartheta_{1,n}-\beta_0) \rightarrow_d J_1^{\textnormal{\textsf{T}}}{\mathbb Z}.$
\item[\textnormal{($b$)}] If $r_n = 3$, let $\beta_0$ and $\beta_0^\dagger$ be the simple and the double root, respectively, where $$G_\beta(\beta_0^\dagger)=0,\ G_{\beta\beta}(\beta_0^\dagger)\eqqcolon 2a,\ \text{ and }\ G_\lambda(\beta_0^\dagger) \eqqcolon B.$$
Then we have for three roots $\vartheta_n$, $\vartheta_n^\pm$, say, the following $$(\sqrt{n}(\vartheta_{n}-\beta_0),n^{1/4}(\vartheta_{n}^\pm-\beta_0^\dagger)) \rightarrow_d (J^{\textnormal{\textsf{T}}}{\mathbb Z},\pm \sqrt{|a^{-1}B^{\textnormal{\textsf{T}}}{\mathbb Z}|}).$$
Moreover, $\sqrt n(\textstyle{\frac{1}{2}}(\vartheta_{n}^+ + \vartheta_{n}^-) - \beta_0^\dagger) = O_p(1)$. If $cB \neq aD$, with $G_{\beta\beta\beta}(\beta_0^\dagger) \eqqcolon 6c$, $G_{\beta\lambda}(\beta_0^\dagger) \eqqcolon D$, then
$$
\sqrt n(\textstyle{\frac{1}{2}}(\vartheta_{n}^+ + \vartheta_{n}^-) - \beta_0^\dagger) \rightarrow_d \frac1{2a}(\frac{c}{a}B-D)^{\textnormal{\textsf{T}}} {\mathbb Z},
$$
and $\sqrt n(\textstyle{\frac{1}{2}}(\vartheta_{n}^+ + \vartheta_{n}^-) - \beta_0^\dagger) = o_p(1)$ otherwise.
\end{itemize}
\end{itemize}
\end{corollary}
Before proceeding, let us briefly comment on the case $r_0 = 2$ where $G(\beta)$ merely touches the horizontal axis. Using an argument similar to \citet[Proof of Thm 1]{dovonon:18}, a second-order Taylor expansion yields, for $\beta$ around the double root $\beta_0^\dagger$,
\[
\frac1{a}G(\beta;\lambda_n) = (\beta - \beta_0^\dagger)^2+{\mathbb Z}_n+ o_p(n^{-1/2}), \quad {\mathbb Z}_n \coloneqq \frac1{a}B^{\textnormal{\textsf{T}}}(\lambda_n-\lambda_0),
\]
with $\sqrt{n}\mathbb{Z}_n \rightarrow_d \frac1{a}B^{\textnormal{\textsf{T}}}{\mathbb Z}$. That is, we have a parabola $\beta \mapsto (\beta-\beta_0^\dagger)^2$ shifted by a random intercept $\mathbb{Z}_n = O_p(n^{-1/2})$. As shown in Figure~\ref{fig:repeated}, a random downward shift (${\mathbb Z}_n < 0$) yields, in addition to the simple root on the left, the dotted line ($r_n = 3$) with two more crossings at $\vartheta_n^\pm$, while (${\mathbb Z}_n > 0$) shifts the parabola above the zero line so that the dashed line with no additional crossing emerges ($r_n = 1$). Because ${\mathbb Z}$ is symmetric, both these events occur asymptotically with equal probability of $\frac1{2}.$ Importantly, although the roots are consistent the number of roots is not. Intuitively, consistency of the repeated roots $\vartheta_n^\pm$ depends on the {\it magnitude} of $\sqrt{|\mathbb{Z}_n|}$, which is $O_p(n^{-1/4})$. On the other hand, the root count $r_n$ is determined by the {\it sign} of $\mathbb{Z}_n$, which, by limiting Gaussianity of $\sqrt{n}\mathbb{Z}_n$, remains a fair coin asymptotically. Put differently, $\lambda \mapsto r_n(\lambda)$ is discontinuous at $\lambda_0$ under multiple roots, so $r(\lambda_n)$ need not to converge to $r_0 = 2$ even though $\vartheta_n^\pm$ are consistent on the $\{r_n = 3\}$ event. An interesting by-product is that the double-roots $\vartheta_n^\pm$ have, with opposite sign, the same leading term of order $n^{-1/4}$. For large $n$, they therefore split symmetrically around the repeated root $\beta_0^\dagger$. Akin to a two-sample Jackknife, averaging thus eliminates this leading {\it bias} term.
\subsection{Inference}\label{sec:inf}
The variance-covariance matrix $\Sigma = \sigma_u^2A^{-1}$, $A = \textnormal{\textsf{E}}[\dot f_t \dot f_t^{\textnormal{\textsf{T}}}]$, of the limiting distribution of $\theta_n$ is not known in closed-form. We therefore resort to numerical derivatives to estimate $A$. For a step size $\ell_n$, the estimator $A_n \coloneqq A_n(\theta_n)$ with $i$,$j$-th element ($1\leq i,j \leq 3$) is given by
\begin{equation}
\begin{split}\label{eq:numA}
[A_{n}(\theta)]_{i,j} \coloneqq \,& \frac1{2 n}(Q_n(\theta+e_i\ell_n+e_j\ell_n)-Q_n(\theta-e_i\ell_n+e_j\ell_n) \\
&\quad\qquad-Q_n(\theta+e_i\ell_n-e_j\ell_n)+Q_n(\theta-e_i\ell_n-e_j\ell_n))/(2\ell_n)^2
\end{split}
\end{equation}
for \(\ell_n \searrow 0\) and $e_i$, $i \in \{1,2,3\}$, denoting some sequence of step-sizes and the $3 \times 1$ unit vector, respectively. Recalling $\hat\sigma_u^2 = \frac1{n} Q_n(\theta_n)$, the following result is obtained:
\begin{corollary}\label{cor:SEs}
If the conditions of Proposition~\textup{\ref{prop:norm}} hold and $\ell_n \rightarrow 0$ as $\sqrt{n}\ell_n\rightarrow \infty$, we have $\Sigma_n \coloneqq \hat\sigma_u^2A_n^{-1} \rightarrow_p \Sigma$, where convergence in probability holds elementwise.
\end{corollary}
The rate imposed on the step-size is common in the literature (see, e.g. \citealp[Theorem 7.4]{newmc:94}, \citealp[Section 2.4]{oh2013simulated} or \citealp[Corollary 2]{mm:25}).
We can use Corollary \ref{cor:SEs} directly to test hypotheses about $\theta_0 \in \Theta$ using Wald-type tests. An important exception is the boundary case $\delta = 0$, where $\gamma$ is a non-identified nuisance parameter; this is similar to \cite{andpol:1994} and \cite{hansen:1996}. As a solution, we propose the following procedure similar to \cite{hansen:17} (see also \citealp{mm:25}): Under the null $H_0$: $\delta = 0$, the ALM reduces to $\pi_t = \psi y_t+u_t$, where $\psi$ can be estimated by the OLS estimator $\tilde\psi_n$, say. Let $\tilde\sigma_n^2 = \frac1{n}\sum_t^n(\pi_t-\tilde\psi_n y_t)^2$ and define $\hat\sigma_n^2 \coloneqq \hat\sigma_n(\gamma_n)$, $\hat\sigma_n(\gamma) \coloneqq \frac1{n}Q_n^\star(\gamma)$, where $Q_n^\star(\cdot)$ is the profiled objective (see Remark~\ref{rem:prof} above). The proposed test statistic is the following {\it sup}{\sf F} statistic:
\begin{align}\label{eq:supF}
{sup}{\sf F}\coloneqq\operatorname*{\textnormal{\textsf{sup}}}\limits_{\gamma \in \Gamma} F_n(\gamma),\quad F_n(\gamma) \coloneqq n\frac{\tilde\sigma_n^2-\hat\sigma_n^2(\gamma)}{\hat\sigma_n^2(\gamma)}.
\end{align}
Note that, by definition of the profiled $\gamma_n$, the statistic takes it maximum at ${\it sup}{\sf F} = n(\tilde\sigma_n-\hat\sigma_n)/\hat\sigma_n$, while, by Frisch-Waugh-Lovell, $F_n(\gamma)$ equals the squared $t$-statistic for the exclusion of $h_{t-1}(\gamma)$ in a regression of $\pi_t$ on $(h_{t-1}(\gamma),y_t)$.
\begin{corollary}\label{cor:sup} Suppose Assumptions \textup{\ref{ass:density}}, \textup{\ref{ass:para}} and \textup{\ref{ass:eight}} hold and $H_0:$ $\delta = 0$, then
\[
{sup}{\sf F} \rightarrow_d \operatorname*{\textnormal{\textsf{sup}}}\limits_{\gamma \in \Gamma} {\mathbb U}(\gamma)^2,
\]
where $ {\mathbb U}(\gamma)$ is a Gaussian process with
$$
\textnormal{\textsf{cov}}[{\mathbb U}(\gamma_1),{\mathbb U}(\gamma_2)] = \frac{\textnormal{\textsf{E}}[h^\star_{t-1}(\gamma_1)h^\star_{t-1}(\gamma_2)]}{\sqrt{\textnormal{\textsf{E}}[h^\star_{t-1}(\gamma_1)^2]\textnormal{\textsf{E}}[h^\star_{t-1}(\gamma_2)^2]}},
$$
where $h^\star_{t}(\cdot)$ is the population residual upon partialling out $y_t$ from $h_t(\cdot)$.
\end{corollary}
Since the limiting process depends on unknown quantities, we follow \cite{hansen:17} and propose the use of a Gaussian multiplier bootstrap to obtain critical values: For $b \in \{1,\dots,B\}$, let $\{\pi_{1,b},\dots,\pi_{n,b}\}$, be a draw of standard normal variates. Use $\pi_{t,b}$ in place of $\pi_t$ to compute the bootstrap equivalent $sup{\sf F}^b$, say, of Eq. \eqref{eq:supF}. For $B$ large enough, we can then approximate the bootstrap $p$-value using $p_{n,B} \coloneqq \frac1{B}\sum_{b=1}^B1\{sup{\sf F} \leq sup{\sf F}^b\}$, and reject at a significance level $\alpha \in (0,1)$ if $p_{n,B} \leq \alpha$. As shown in the appendix $p_{n,B}$ is consistent for the true $p$-value as $n$ and $B$ diverge.
Another interesting question is how to conduct inference with respect to the equilibria $\beta$. Here, we propose uniform confidence bands for the unknown curve $\beta \mapsto G(\beta)$. To fix ideas, note that, by the functional Delta method, we have in $\ell^\infty([0,1])$, $\sqrt{n}(G(\beta,\lambda_n) - G(\beta,\lambda_0)) \rightsquigarrow {\mathbb G}(\beta),$
where ${\mathbb G}(\beta)$ is a Gaussian process with kernel $\textnormal{\textsf{cov}}[{\mathbb G}(\beta_1),{\mathbb G}(\beta_2)] = F_\lambda(\beta_1)^{\textnormal{\textsf{T}}}\Omega F_\lambda(\beta_2)$. We can therefore compute uniform $1-\alpha$ confidence intervals as summarized by the following Corollary:
\begin{corollary}\label{cor:CI} For a consistent estimator $s_n^2(\beta)$ of $s^2(\beta) \coloneqq F_\lambda(\beta)^{\textnormal{\textsf{T}}}\Omega F_\lambda(\beta)$ such that $\operatorname*{\textnormal{\textsf{sup}}}_\beta|s_n(\beta)-s(\beta)| = o_p(1)$, fix a significance level $\alpha \in (0,1)$ and suppose Assumptions \textup{\ref{ass:density}}, \textup{\ref{ass:para2}}, and \textup{\ref{ass:eight}} hold, then $$\mathbb{P}\{\forall \beta \in [0,1]: G(\beta) \in [G_n(\beta) \pm c_\alpha n^{-1/2}s_n(\beta)]\} \rightarrow 1-\alpha,$$ where $c_\alpha$ is such that $\mathbb{P}\{\operatorname*{\textnormal{\textsf{sup}}}\limits_{\beta \in [0,1]} \frac{|{\mathbb G}(\beta)|}{s(\beta)} \leq c_\alpha\} = 1-\alpha$.
\end{corollary}
If preferred, a non-studentized, {\it percentile} version of the confidence bands could be obtained by letting $[G_n(\beta) \pm \tilde c_\alpha n^{-1/2}]$, with $\tilde c_{\alpha}$ such that $\mathbb{P}\{\operatorname*{\textnormal{\textsf{sup}}}_{\beta} |{\mathbb G}(\beta)| \leq \tilde c_\alpha\} = 1-\alpha$. Note that, contrary to pointwise confidence intervals at selected roots $\beta$, our confidence bands are constructed for the function $\beta \mapsto G_n(\beta;\lambda_0)$ itself, exploiting smoothness in $\lambda$, and are thus unaffected by potentially repeated roots where $G_\beta$ vanishes.
Since the limiting distribution is unknown and non-pivotal, to obtain $c_\alpha$, we suggest to use a multiplier bootstrap based on the first order approximation:
\[
\sqrt{n}(G_n(\beta)-G(\beta)) = \frac1{\sqrt n}\sum_{t=1}^n m_t(\beta;\kappa_0)+o_p(1),
\]
where
\[
m_t(\beta;\kappa) \coloneqq -F_\lambda(\beta;\lambda)^{\textnormal{\textsf{T}}}\phi_t(\kappa), \quad \kappa \coloneqq (\lambda^{\textnormal{\textsf{T}}},\gamma)^{\textnormal{\textsf{T}}},
\]
and $\frac1{\sqrt n}\sum_t \phi_t(\kappa_0)$ is the influence function of $\sqrt{n}(\lambda_n-\lambda_0)$. Thus, for a draw $b \in \{1,\dots,B\}$ of standard Gaussian multipliers $\{\xi_{1,b},\dots,\xi_{n,b}\}$, and, conditionally on the data, we have
\[
\frac1{\sqrt n}\sum_{t=1}^n \xi_{t,b} m_t(\beta;\kappa_0) \rightsquigarrow_{\mathbb P} {\mathbb G}(\beta) \quad \textup{in} \quad \ell^\infty([0,1]),
\]
where $ \rightsquigarrow_{\mathbb P}$ means that weak convergence holds in probability with respect to the data-generating measure. In practice, neither $\kappa_0$ nor $\phi_t(\kappa_0)$ is observed. As shown in Appendix~\ref{app:inf}, we can however construct $\phi_{t,n}(\kappa_n)$ such that
$\frac{1}{n}\sum_{t=1}^n\|\phi_{t,n}(\kappa_n)-\phi_t(\kappa_0)\|^2=o_p(1)$, and define the feasible linearization
$m_{t,n}(\beta)\coloneqq -F_\lambda(\beta;\lambda_n)^\top\phi_{t,n}(\kappa_n)$. We then simulate the critical value from the studentized bootstrap statistic
\[
T_{n,b} \coloneqq \operatorname*{\textnormal{\textsf{sup}}}_{\beta\in[0,1]}\frac{\big|\frac{1}{\sqrt n}\sum_{t=1}^n \xi_{t,b} m_{t,n}(\beta)\big|}{s_n(\beta)},
\qquad
s_n(\beta)^2\coloneqq \frac{1}{n}\sum_{t=1}^n m_{t,n}(\beta)^2,
\]
and take $c_{\alpha,n,B}$ as the empirical $1-\alpha$ quantile of $\{T_{n,b}\}_{b=1}^B$. As verified in the appendix, if $B \rightarrow \infty$ and $n \rightarrow \infty$, $c_{\alpha,n,B}\rightarrow_p c_\alpha$, yielding the uniform bands in Corollary~\ref{cor:CI}.
\section{Monte Carlo simulation}\label{sec:mc}
In this section, we present finite-sample evidence based on simulated data. We consider two scenarios: Scenario~(A) with three behavioural equilibria ($r_0 = 3$) and Scenario~(B) with two equilibria ($r_0 = 2$). The Monte Carlo dgp is defined by Eqs.~\eqref{eq:AR1}, \eqref{eq:learning}, and \eqref{eq:ALM}, with {\sf IID} Gaussian innovations $v_t = (u_t,\varepsilon_t)^{\textnormal{\textsf{T}}} \sim \mathcal N(0,\textup{\sf diag}(\sigma_u^2,\sigma_\varepsilon^2))$. For Scenario~(A), the parameter values are inspired by our empirical analysis. In particular, the parameter vector estimated by NLS is set to $\theta = (0.076,0.998,0.090)^{\textnormal{\textsf{T}}}$, the AR(1) dynamics~\eqref{eq:AR1} are parametrised via $(a,\rho) = (-0.02,0.93)$, and, for the innovations, we set $(\sigma_u,\sigma_\varepsilon) = (0.44,0.76)$. This configuration gives rise to three behavioural equilibria $\beta = (0.174,0.856,0.999)^{\textnormal{\textsf{T}}}$. The implied values of $\lambda = (\delta,\psi,\rho,\sigma_u^2,\sigma_\varepsilon^2)^{\textnormal{\textsf{T}}}$ are broadly in line with those in \citet[Sec~4.2]{hz14}. For Scenario~(B), we keep the value of the gain and the AR(1) dynamics, and choose
$\lambda = (0.95,0.2734,0.9,1,1)^{\textnormal{\textsf{T}}}$, which gives rise to a simple root at $\beta = 0.9766$ and a repeated root at $\beta^\dagger = 0.5551$.
\begin{figure}[!h]
\centering
\begin{minipage}[c]{.7\textwidth}
\begin{tikzpicture}
\begin{axis}[
width=\linewidth,
height=5cm,
xlabel={$\beta$},
ylabel={$G(\beta)$},
legend pos=south east,
legend cell align=right,
legend style={draw=none, at={(rel axis cs:0.9,.05)}, fill=none,font=\footnotesize},
table/col sep=space,
grid style={line width=.1pt, draw=gray!20},
major grid style={line width=.2pt, draw=gray!40},
scaled y ticks=false,
yticklabel style={/pgf/number format/fixed, /pgf/number format/precision=4},
enlarge x limits=false,
xticklabels={},ytick={0},
yticklabels={0},
xlabel near ticks,
ylabel near ticks,
label style={font=\footnotesize},
tick label style={font=\scriptsize},
xlabel style={yshift=2pt},
ylabel style={xshift=2pt},
xmin = 0.1, xmax = 1.01,
ymin = -.11, ymax = .13,
]
\draw[red, line width=.5pt]
(axis cs:\pgfkeysvalueof{/pgfplots/xmin}, 0)
-- (axis cs:\pgfkeysvalueof{/pgfplots/xmax}, 0);
Thicker grey dotted vertical lines
\draw[gray, line width=.5pt]
(axis cs:0.9986868, \pgfkeysvalueof{/pgfplots/ymin})
-- (axis cs:0.9986868, \pgfkeysvalueof{/pgfplots/ymax});
\draw[gray, line width=.5pt]
(axis cs:0.9206429, \pgfkeysvalueof{/pgfplots/ymin})
-- (axis cs:0.9206429, \pgfkeysvalueof{/pgfplots/ymax});
\draw[gray, line width=.5pt]
(axis cs:0.1345104, \pgfkeysvalueof{/pgfplots/ymin})
-- (axis cs:0.1345104, \pgfkeysvalueof{/pgfplots/ymax});
\draw[gray, dashed, line width=.5pt]
(axis cs:0.5551, \pgfkeysvalueof{/pgfplots/ymin})
-- (axis cs:0.5551, \pgfkeysvalueof{/pgfplots/ymax});
\draw[gray, dashed,line width=.5pt]
(axis cs:.9766 , \pgfkeysvalueof{/pgfplots/ymin})
-- (axis cs:.9766 , \pgfkeysvalueof{/pgfplots/ymax});
\addplot+[smooth, no marks, no marks, thick, black] table[x=betas, y=G1] {gfun.txt}; \addlegendentry{(A)}
\addplot+[smooth, dashed, no marks, black] table[x=betas, y=G2] {gfun.txt}; \addlegendentry{(B)}
\end{axis}
\end{tikzpicture}
\end{minipage}
\begin{minipage}[c]{.275\textwidth}\vspace*{-.5cm}
\caption{The solid and dashed lines indicate the true $G(\beta)$ for scenario (A) and (B), respectively. The $r_0 = 3$ and $r_0=2$ roots are indicate with grey vertical lines.}\label{fig:mcfun}
\end{minipage}
\end{figure}
\begin{table}[!hp]
\centering
\scriptsize
\setlength{\tabcolsep}{3pt}
\resizebox{\linewidth}{!}{
\begin{tabular}{lrcccrrccc}
\toprule
(A) & & $\gamma = 0.076$ & $\delta = 0.998$ & $\psi = 0.090$ & & & $\beta_1 = 0.174$ & $\beta_2 = 0.856$ & $\beta_3 = 0.999$ \\
\midrule
250 & mean & 0.0754 & 0.9696 & 0.0929 & & $\mathbb{P}\{r=1\}$ & & {\bf 0.4085} & \\
& bias & -0.0006 & -0.0285 & 0.0029 & & mean & & 0.2250 & \\
& sd & 0.0171 & 0.0551 & 0.0160 & & sd & & 0.2911 & \\[-2pt]
\cmidrule(r){8-10} \\[-8.4pt]
& size & 0.0595 & 0.0415 & 0.0525 & & $\mathbb{P}\{r=3\}$ & & {\bf 0.5915 } & \\
& cover & 0.9335 & & & & mean & 0.1924 & 0.8516 & 0.9961 \\
& cover-s & 0.9085 & & & & sd & 0.0832 & 0.0929 & 0.0086 \\[-2pt]
\cmidrule(r){1-5} \cmidrule(r){7-10}
500 & mean & 0.0760 & 0.9860 & 0.0911 & & $\mathbb{P}\{r=1\}$ & & {\bf 0.2310} & \\
& bias & 0.0000 & -0.0121 & 0.0011 & & mean & & 0.1967 & \\
& sd & 0.0123 & 0.0240 & 0.0110 & & sd & & 0.2434 & \\[-2pt]
\cmidrule(r){8-10} \\[-8.4pt]
& size & 0.0565 & 0.0440 & 0.0500 & & $\mathbb{P}\{r=3\}$ & & {\bf 0.7690} & \\
& cover & 0.9290 & & & & mean & 0.1863 & 0.8547 & 0.9963 \\
& cover-s & 0.9280 & & & & sd & 0.0692 & 0.0802 & 0.0072 \\[-2pt]
\cmidrule(r){1-5} \cmidrule(r){7-10}
1,000 & mean & 0.0764 & 0.9935 & 0.0905 & & $\mathbb{P}\{r=1\}$ & & {\bf 0.0695} & \\
& bias & 0.0004 & -0.0046 & 0.0005 & & mean & & 0.1229 & \\
& sd & 0.0086 & 0.0104 & 0.0077 & & sd & & 0.0292 & \\[-2pt]
\cmidrule(r){8-10} \\[-8.4pt]
& size & 0.0455 & 0.0285 & 0.0495 & & $\mathbb{P}\{r=3\}$ & & {\bf 0.9305} & \\
& cover & 0.9350 & & & & mean & 0.1778 & 0.8598 & 0.9972 \\
& cover-s & 0.9490 & & & & sd & 0.0494 & 0.0600 & 0.0050 \\[-2pt]
\cmidrule(r){1-5} \cmidrule(r){7-10}
2,000 & mean & 0.0758 & 0.9962 & 0.0904 & & $\mathbb{P}\{r=1\}$ & & {\bf 0.0120} & \\
& bias & -0.0002 & -0.0018 & 0.0004 & & mean & & 0.1298 & \\
& sd & 0.0060 & 0.0048 & 0.0054 & & sd & & 0.0224 & \\[-2pt]
\cmidrule(r){8-10}
& size & 0.0520 & 0.0345 & 0.0480 & & $\mathbb{P}\{r=3\}$ & & \textbf{0.9880} & \\
& cover & 0.9330 & & & & mean & 0.1759 & 0.8582 & 0.9983 \\
& cover-s & 0.9505 & & & & sd & 0.0357 & 0.0441 & 0.0030 \\
\midrule
(B) & & $\gamma = 0.076$ & $\delta = 0.95 $ & $\psi = 0.2734 $ & & & $\beta^\dagger$ & $\beta^\dagger=0.55$ & $\beta=0.977$ \\
\midrule
250 & mean & 0.0759 & 0.9361 & 0.2772 & & $\mathbb{P}\{r=1\}$ & & {\bf 0.7550} & \\
& bias & -0.0001 & -0.0139 & 0.0039 & & mean & & 0.6382 & \\
& sd & 0.0156 & 0.0507 & 0.0310 & & sd & & 0.3632 & \\[-2pt]
\cmidrule(r){8-10} \\[-8.4pt]
& size & 0.0540 & 0.0450 & 0.0460 & & $\mathbb{P}\{r=3\}$ & & {\bf 0.2450} & \\
& cover & 0.9270 & & & & mean & 0.3429 & 0.7671 & 0.9755 \\
& cover-s & 0.9125 & & & & sd & 0.0954 & 0.0906 & 0.0242 \\
& & & & & & & \multicolumn{2}{c}{$\textstyle{\frac{1}{2}}(\beta_1+\beta_2)$} & \\[-1pt] \cmidrule(r){8-9}
& & & & & & mean & \multicolumn{2}{c}{0.5550} & \\
& & & & & & sd & \multicolumn{2}{c}{0.0255} & \\[-2pt]
\cmidrule(r){1-5} \cmidrule(r){7-10}
500 & mean & 0.0762 & 0.9438 & 0.2746 & & $\mathbb{P}\{r=1\}$ & & {\bf0.6525} & \\
& bias & 0.0002 & -0.0062 & 0.0012 & & mean & & 0.7166 & \\
& sd & 0.0113 & 0.0261 & 0.0216 & & sd & & 0.3428 & \\[-2pt]
\cmidrule(r){8-10}
& size & 0.0520 & 0.0490 & 0.0515 & & $\mathbb{P}\{r=3\}$ & & {\bf0.3475} & \\
& cover & 0.9255 & & & & mean & 0.3764 & 0.7572 & 0.9641 \\
& cover-s & 0.9205 & & & & sd & 0.0771 & 0.0804 & 0.0239 \\
& & & & & & & \multicolumn{2}{c}{$\textstyle{\frac{1}{2}}(\beta_1+\beta_2)$} & \\[-1pt] \cmidrule(r){8-9}
& & & & & & mean & \multicolumn{2}{c}{0.5668} & \\
& & & & & & sd & \multicolumn{2}{c}{0.0236} & \\[-2pt]
\cmidrule(r){1-5} \cmidrule(r){7-10}
1,000 & mean & 0.0766 & 0.9477 & 0.2737 & & $\mathbb{P}\{r=1\}$ & & {\bf0.5615} & \\
& bias & 0.0006 & -0.0023 & 0.0004 & & mean & & 0.8238 & \\
& sd & 0.0077 & 0.0138 & 0.0152 & & sd & & 0.2913 & \\[-2pt]
\cmidrule(r){8-10}
& size & 0.0470 & 0.0400 & 0.0445 & & $\mathbb{P}\{r=3\}$ & & {\bf0.4385} & \\
& cover & 0.9385 & & & & mean & 0.3893 & 0.7468 & 0.9637 \\
& cover-s & 0.9430 & & & & sd & 0.0680 & 0.0793 & 0.0168 \\
& & & & & & & \multicolumn{2}{c}{$\textstyle{\frac{1}{2}}(\beta_1+\beta_2)$} & \\[-1pt] \cmidrule(r){8-9}
& & & & & & mean & \multicolumn{2}{c}{0.5680} & \\
& & & & & & sd & \multicolumn{2}{c}{0.0166} & \\[-2pt]
\cmidrule(r){1-5} \cmidrule(r){7-10}
2,000 & mean & 0.0760 & 0.9487 & 0.2739 & & $\mathbb{P}\{r=1\}$ & & {\bf0.5015} & \\
& bias & 0.0000 & -0.0013 & 0.0006 & & mean & & 0.9253 & \\
& sd & 0.0054 & 0.0087 & 0.0105 & & sd & & 0.1863 & \\[-2pt]
\cmidrule(r){8-10}
& size & 0.0585 & 0.0445 & 0.0480 & & $\mathbb{P}\{r=3\}$ & & {\bf0.4985} & \\
& cover & 0.9415 & & & & mean & 0.4066 & 0.7263 & 0.9659 \\
& cover-s & 0.9425 & & & & sd & 0.0608 & 0.0752 & 0.0136 \\
& & & & & & & \multicolumn{2}{c}{$\textstyle{\frac{1}{2}}(\beta_1+\beta_2)$} & \\[-1pt] \cmidrule(r){8-9}
& & & & & & mean & \multicolumn{2}{c}{0.5664} & \\
& & & & & & sd & \multicolumn{2}{c}{0.0133} & \\
\bottomrule
\end{tabular}
}\caption{Simulation results for scenario (A) in the upper panel with $r_0 =3$ roots and for scenario (B) in the lower panel with $r_0=2$ roots based on 2,000 Monte Carlo repetitions. }\label{tab:MC}
\end{table}
We report, for the NLS estimator $\theta_n$ of $\theta = (\gamma,\delta,\psi)^{{\textnormal{\textsf{T}}}}$, the mean, bias, standard deviation, and the empirical size of two-sided $t$-tests (using, in accordance with Corollary~\ref{cor:SEs}, numerical standard errors) at a nominal 5\% level. We also present the empirical coverage of the percentile ({\it cover}) and studentized ({\it cover-s}) confidence bands (see Corollary~\ref{cor:CI}) for the curves $\beta \mapsto G(\beta)$ plotted in Figure~\ref{fig:mcfun}, based on $B = 4{,}999$ bootstrap replications at a 95\% coverage level. Finally, we report the mean and standard deviation of the estimated equilibria, together with the observed frequency at which the true number of roots is recovered, i.e.\ $\mathbb{P}\{r = r_0\}$.
The findings are summarized in Table~\ref{tab:MC}. In both scenarios, the NLS estimator appears consistent and the normal approximation for the $t$-statistics is fairly accurate, while the empirical coverage of the studentized confidence bands is slightly below the nominal 95\% level for smaller $n$. Moreover, in line with Corollary~\ref{cor:roots}, we observe that the estimated number of roots $r_n$ either converges to the true value $r_0$ or, in Scenario~(B), oscillates between $r = 1$ and $r = 3$ with approximately equal probability. In the latter case, the convergence of the estimates associated with the repeated root appears to be an order of magnitude slower than the $\sqrt{n}$-rate observed for simple roots, while averaging the two nearby estimates corresponding to the repeated root improves performance.
\section{Empirical application}\label{sec:emp}
To illustrate the empirical relevance of our theoretical results, we estimate the NKPC given by Eqs.~\eqref{eq:AR1}, \eqref{eq:learning}, and \eqref{eq:ALM} using quarterly US\ data from 1960:Q1 to 2019:Q4 ($n=240$). Inflation $\pi_t$ is measured by CPI inflation. Since marginal costs are proportional to the output gap or the real labour share (see, e.g., \citealp{Woodford2003} or \citealp[Section 2.1]{MavroeidisPlagborgMollerStock2014}), we consider two standard slack measures: ({\sf A}) the output gap measure from the Congressional Budget Office and ({\sf B}) the percentage change in real unit labour costs (see also \citealp{gali1999inflation} or \citealp{sbordone2002prices}). The three series are plotted in Figure~\ref{fig:empraw}.
\begin{figure}[!h]
\centering
\begin{tikzpicture}
\begin{axis}[
name=bottom,
width=13cm,
height=5cm,
at={(0,0)},
anchor=south west,
xmin=1960,
xmax=2020,
axis lines=box,
xtick={1960,1970,1980,1990,2000,2010,2020},
xticklabel style={
/pgf/number format/fixed,
/pgf/number format/precision=0,
/pgf/number format/set thousands separator={}
},
xlabel={},
ylabel={{\it output gap}},
y label style={color=blue},
yticklabel style={color=blue},
grid=none,
table/col sep=space
]
\addplot+[no marks, thick, blue] table[
x expr=1960 + \coordindex/4,
y=y
] {raw.txt};
\end{axis}
\begin{axis}[
name=top,
width=13cm,
height=4cm,
at={(bottom.north west)},
anchor=south west,
xmin=1960,
xmax=2020,
axis lines=box,
xtick={1960,1970,1980,1990,2000,2010,2020},
xticklabels={},
xlabel={},
ylabel={$\pi$},
grid=none,
table/col sep=space
]
\addplot+[no marks, black, thick] table[
x expr=1960 + \coordindex/4,
y=pi
] {raw.txt};
\end{axis}
\begin{axis}[
width=13cm,
height=4cm,
at={(bottom.south west)},
anchor=south west,
xmin=1960,
xmax=2020,
axis y line*=right,
axis x line=none,
ylabel={{\it unit labour costs}},
y label style={color=pink},
yticklabel style={color=pink},
ytick pos=right,
grid=none,
table/col sep=space
]
\addplot+[no marks, thick, pink] table[
x expr=1960 + \coordindex/4,
y=ya
] {raw.txt};
\end{axis}
\end{tikzpicture}
\caption{Percentage change of the CPI (upper panel) and, in the lower panel, the output gap and the percentage change of unit labour costs .}\label{fig:empraw}
\end{figure}
We employ a preliminary step to obtain starting values: In particular, we use an initial grid search over $\gamma \in (0,0.1)$ in order to obtain starting values, treating the initial value $\operatorname*{\mathfrak{a}}$ of the learning recursion as a parameter to be estimated. NLS estimation is carried out using the {\sf R} {\sf optim}-routine with default settings, restricting the parameter space for $\delta$ to $(0,1)$. The upper panel of Table~\ref{tab:emp} reports the NLS estimation results, while the lower panel summarizes OLS estimation of the AR(1) models for the two different proxies of the driving variable; 95\% confidence intervals based on numerical derivatives are given beneath. The NLS estimates of specification {\sf A} (output gap) and {\sf B} (unit labour costs) are in line with the empirical macro literature (see, e.g., \citealp{milani:07}, \citealp{chev:10}, \citealp{Lansing2009}, or \citealp{HommesMavromatisOzdenZhu2023}). The estimated learning gains imply that agents place most weight on roughly the last 3 to 5 years of data for specifications {\sf A} and {\sf B}, respectively. The estimate of the discount factor on expected inflation (viz. $\delta$) is large. This confirms that forward-looking expectations play a quantitatively important role in the Phillips curve and that inflation dynamics are highly persistent. The slope parameter is positive and statistically significant, but small in magnitude, consistent with the `flattened' NKPC found in much of the empirical literature.
\begin{table}[ht]
\centering
\begin{tabular}{lcc}
\toprule
& {\sf A} & {\sf B} \\
\midrule
$\gamma$ & 0.07602311 & 0.04630777 \\[-5pt]
& {\scriptsize [0.07292689,\ 0.07911934]} & {\scriptsize [0.04489319,\ 0.04772236]}\\
$\delta$ & 0.998593 & 0.90173298 \\[-5pt]
& {\scriptsize [0.9981619,\ 0.9990241]} & {\scriptsize [0.89704627,\ 0.9064197]}\\
$\psi$ & 0.08897832 & 0.12586481\\[-5pt]
& {\scriptsize [0.0873151,\ 0.09064154]} & {\scriptsize [0.12197922,\ 0.12975041]}\\
$\sigma_u$ &0.44388612& 0.4673309 \\[-1pt] \cmidrule(r){2-3}
$a$ & -0.0259546 & 0.5660696 \\[-5pt]
& {\scriptsize [-0.1253596,\ 0.0734504]} & {\scriptsize [0.4009362,\ 0.731203]}\\
$\rho$ &0.9375791& 0.15752518\\[-5pt]
& {\scriptsize [0.8929128,\ 0.9822454]} & {\scriptsize [0.03131037,\ 0.28374]}\\
$\sigma_\varepsilon$&0.7613624& 1.1100973\\
\bottomrule
\end{tabular}
\caption{Point estimates and 95\% confidence intervals based on numerical standard errors.}\label{tab:emp}
\end{table}
\begin{figure}[!h]
\centering
\begin{tikzpicture}
\begin{axis}[
width=13cm,
height=7cm,
xmin=0, xmax=1,
axis lines=box,
xlabel={$\beta$},
ylabel={},
grid=none,
table/col sep=tab,
legend style={
draw=none,
fill=none,
at={(0.97,0.03)},
anchor=south east
},
legend cell align=left
]
\addplot[gray, dashed, forget plot]
coordinates {(0,0) (1,0)};
\addplot[name path=G1up, draw=none, forget plot]
table[x=grid, y=G1up]{ci.txt};
\addplot[name path=G1low, draw=none, forget plot]
table[x=grid, y=G1low]{ci.txt};
\addplot[
fill=pink,
fill opacity=0.25,
draw=none,
forget plot
]
fill between[of=G1low and G1up];
\addplot[name path=G2up, draw=none, forget plot]
table[x=grid, y=G2up]{ci.txt};
\addplot[name path=G2low, draw=none, forget plot]
table[x=grid, y=G2low]{ci.txt};
\addplot[
fill=blue,
fill opacity=0.1,
draw=none,
forget plot
]
fill between[of=G2low and G2up];
\addplot[pink, thick]
table[x=grid, y=G1]{ci.txt};
\addplot[blue, thick]
table[x=grid, y=G2]{ci.txt};
\legend{{\it unit labour costs}, {\it output gap}}
\end{axis}
\end{tikzpicture}
\caption{Plot of $\beta \mapsto G_n(\beta)$ together with 90\% confidence bands based on $B=$ 4,999 bootstrap replications.}\label{fig:empci}
\end{figure}
In specification {\sf A}, the estimated parameters imply three behavioural equilibria $\beta_1 = 0.1914$, $\beta_2 = 0.8316$, $\beta_3= 0.9995$; in specification {\sf B} we find a unique equilibrium $\beta = 0.01341$. The corresponding $\beta \mapsto G_n(\beta)$ functions are shown in Figure~\ref{fig:empci}. The key difference is the persistence of the driving variable: $\rho$ is close to one in {\sf A} but relatively small in {\sf B}. This pattern is consistent with the theoretical properties of the map $\beta \mapsto G(\beta;\lambda)$, which shows that, holding other parameters fixed, $\rho \rightarrow 0$ leads to a single equilibrium with $\beta \rightarrow 0$. The low persistence of unit labour costs is sufficient to support only a single, less persistent equilibrium, and the feedback between expectations and realized inflation is too weak to generate a high-persistence belief regime. Finally, in case of specification {\sf A}, $\beta_2$ is under the criterion of \citet[Prop 4]{hz14} not stable; i.e. an agent that uses asymptotically all data ($\gamma \rightarrow 0$) is not able to learn $\beta_2$ in the limit. The empirical model thus points to two stable belief regimes for inflation persistence: a low and a high persistence equilibrium.
\section{Conclusion}\label{sec:conc}
We have developed estimation and inference methods for a New Keynesian Phillips
curve with constant-gain learning and potentially multiple behavioural
equilibria. Under mild conditions the resulting model is geometrically
ergodic, and a nonlinear least squares estimator of the structural parameters is
strongly consistent and asymptotically normal. We explain how to conduct inference with respect to structural parameters, including the equilibria, and apply these techniques to U.S.\ inflation data, revealing multiple belief regimes in the data. This paper is only a starting point and several extensions are left for future work: A natural next step is to develop tools for
multivariate versions of the model, along the lines of
\citet{milani:07} or \citet{HommesMavromatisOzdenZhu2023}, in order to allow for
richer joint dynamics of macroeconomic variables (see also the discussion in \citealp{MavroeidisPlagborgMollerStock2014}). It would also be useful to
relax some of our regularity conditions, to exploit panel data to increase
effective sample size, and to adapt the analysis to decreasing-gain learning
rules using results from earlier work on the econometrics of models with adaptive learning like \citealp{chev:10}, \cite{chev:17}, \cite{chrismass:18}, \cite{mayer:22}, or \cite{mm:25}.
\addcontentsline{toc}{section}{References}
\bibliography{bibl}
\newpage\clearpage