# Semiparametric Functional Factor Models with
Bayesian Rank Selection

Daniel R. Kowal and Antonio Canale Dobelman Family Assistant Professor,
Department of Statistics, Rice University (daniel.kowal@rice.edu).Associate
Professor, Department of Statistics, University of Padova
(antonio.canale@unipd.it).

###### Abstract

Functional data are frequently accompanied by a parametric template that
describes the typical shapes of the functions. However, these parametric
templates can incur significant bias, which undermines both utility and
interpretability. To correct for model misspecification, we augment the
parametric template with an infinite-dimensional nonparametric functional
basis. The nonparametric basis functions are learned from the data and
constrained to be orthogonal to the parametric template, which preserves
distinctness between the parametric and nonparametric terms. This distinctness
is essential to prevent functional confounding, which otherwise induces severe
bias for the parametric terms. The nonparametric factors are regularized with
an ordered spike-and-slab prior that provides consistent rank selection and
satisfies several appealing theoretical properties. The versatility of the
proposed approach is illustrated through applications to synthetic data, human
motor control data, and dynamic yield curve data. Relative to parametric and
semiparametric alternatives, the proposed semiparametric functional factor
model eliminates bias, reduces excessive posterior and predictive uncertainty,
and provides reliable inference on the effective number of nonparametric
terms—all with minimal additional computational costs.

KEYWORDS:  factor analysis; nonparametric regression; shrinkage prior; spike-
and-slab prior; yield curve

##  1 Introduction

###  1.1 Setting and goals

As high-resolution monitoring and measurement systems generate vast quantities
of complex and highly correlated data, _functional data analysis_ has become
increasingly vital for many scientific, medical, business, and industrial
applications. Functional data are (noisy) realizations of random functions
$\\{Y_{i}\\}_{i=1}^{n}$ observed over a continuous domain, such as time,
space, or wavelength, and exhibit a broad variety of shapes. The concurrence
of _complex_ and _voluminous_ data prompts the common use of nonparametric
models for functional data. Yet in many cases, there is valuable information
regarding the functional form of $Y_{i}$. Template curves that describe the
shape of $Y_{i}$ are derived from fundamental scientific laws or motivated by
extensive empirical studies, and often are the focal point of the analysis.
Prominent examples include human motor control (Ramsay,, 2000; Goldsmith and
Kitago,, 2016), interest rates (Nelson and Siegel,, 1987; Diebold and Li,,
2006), and basal body temperature (Scarpa and Dunson,, 2009, 2014; Canale et
al.,, 2017).

Our goal is to construct a functional data modeling framework that
simultaneously (i) incorporates parametric templates in a coherent and
interpretable manner, (ii) maintains the modeling flexibility of nonparametric
methods, and (iii) provides computationally scalable inference with reliable
uncertainty quantification. The approach is fully Bayesian, accompanied by a
highly efficient Markov chain Monte Carlo (MCMC) algorithm for posterior and
predictive inference, and equally applicable to both densely-observed and
sparsely- or irregularly-sampled functional data.

The parametric templates are represented as a spanning set
$\mathcal{H}_{0}=\mbox{span}\\{g_{1}(\cdot;\gamma),\ldots,g_{L}(\cdot;\gamma)\\}$
of functions $\\{g_{\ell}(\cdot;\gamma)\\}_{\ell=1}^{L}$ known up to $\gamma$.
Any function belonging to $\mathcal{H}_{0}$ is a linear combination of
$\\{g_{\ell}\\}$; the corresponding coefficients and the parameters $\gamma$
must be learned. Important examples are presented in Table 1. The linear basis
is routinely used for longitudinal data analysis and here is equivalent to a
random slope model (Molenberghs and Verbeke,, 2000). A change in slope,
$(\tau-\gamma)_{+}$ with $(x)_{+}=\max\\{0,x\\}$, is useful for modeling
structural changes, such as a change in disease transmissions due to policy
interventions (Wagner et al.,, 2020). Cosinor functions model circadian
rhythms (Mikulich et al.,, 2003) and other periodic behaviors (Welham et al.,,
2006). Biphasic curves offer utility in modeling basal body temperature of
women during the menstrual cycle (Scarpa and Dunson,, 2009, 2014; Canale et
al.,, 2017). In general, interest centers on learning the linear coefficients
associated with each $g_{\ell}$, the nonlinear parameters $\gamma$, and an
adequate yet interpretable model for the functions $Y_{i}$.

Linear | Linear change | Cosinor | Biphasic
---|---|---|---
$\\{1,\tau\\}$ | $\\{1,\tau,(\tau-\gamma)_{+}\\}$ | $\\{1,\sin(2\pi\tau/\gamma),\cos(2\pi\tau/\gamma)\\}$ | $\\{1,\exp(\gamma\tau)/\\{1+\exp(\gamma\tau)\\}\\}$
Table 1: Examples of parametric templates, or spanning sets, for
$\mathcal{H}_{0}$.

The advantages of the parametric templates are clear: they incorporate domain
knowledge, lend interpretability to the model, and often produce low-variance
estimators relative to nonparametric alternatives. These templates are
restrictive by design, and therefore can incur significant bias and other
model misspecifications. Such effects erode model interpretability and can
induce variance inflation. For illustration, we present model fits for two
datasets in Figure 1 using a parametric model and our proposed semiparametric
alternative; details and analyses of these data are in Section 6. In both
instances, the templates capture the general shape of the data. However,
regions of substantial bias are present, which produce uniformly inflated
prediction bands over the domain. By comparison, the proposed approach
preserves the essential shape of the curves, yet corrects the bias and shrinks
the prediction bands appropriately—and crucially does so without overfitting.

Figure 1: Posterior expectations and 95% simultaneous prediction bands for the
parametric (PFFM, top) and semiparametric functional factor model (SFFM,
bottom) for one replicate of the pinch force data (left; see Section 6.1) and
one yield curve (right; see Section 6.2). The proposed SFFM corrects the bias
of the PFFM and offers narrower prediction bands.

###  1.2 Overview of the proposed approach

The _semiparametric functional factor model_ (SFFM) bridges the gap between
parametric and nonparametric functional data models. The SFFM augments a
parametric template $\mathcal{H}_{0}$ with a nonparametric and infinite-
dimensional basis expansion for the functions $Y_{i}\in L^{2}(\mathcal{T})$:

$Y_{i}(\tau)=\sum_{\ell=1}^{L}\alpha_{\ell,i}g_{\ell}(\tau;\gamma)+\sum_{k=1}^{\infty}\beta_{k,i}f_{k}(\tau),\quad\tau\in\mathcal{T},$
(1)

where $\\{f_{k}\\}$ are unknown functions, $\\{\alpha_{\ell,i}\\}$ and
$\\{\beta_{k,i}\\}$ are unknown factors for the parametric and nonparametric
components, respectively, and $\mathcal{T}\subset\mathbb{R}^{d}$ is the
domain, usually with $d=1$ for curves or $d=2$ for images. Without the
nonparametric terms, (1) is a _parametric functional factor model_ (PFFM),
which serves as our parametric baseline.

Our implementation of (1) is fully Bayesian with three unique and essential
features.

First, the nonparametric basis $\\{f_{k}\\}$ is treated as unknown. The SFFM
pools information across all functions $\\{Y_{i}\\}$ to learn the key
functional features in the data—in particular, the systemic biases unresolved
by $\mathcal{H}_{0}$. The functions $\\{f_{k}\\}$ are learned jointly with the
parametric terms $\gamma$ and $\\{\alpha_{\ell,i}\\}$ and the nonparametric
factors $\\{\beta_{k,i}\\}$, so uncertainty about $\\{f_{k}\\}$ is
automatically absorbed into the joint posterior distribution.

Second, the nonparametric basis $\\{f_{k}\\}$ is constrained to be orthogonal
to the templates $\\{g_{\ell}\\}$. This important constraint enforces
distinctness between the parametric and nonparametric components of the SFFM.
We show that without such a constraint, inference on the parametric terms
becomes severely biased and loses interpretability due to _functional
confounding_. The nonredundancy of $\\{f_{k}\\}$ is also essential for valid
inference on the effective number of nonparametric terms. Our chosen
constraint produces computational simplifications that improve both the
algorithmic efficiency and the ease of implementation of the MCMC algorithm.

Third, the nonparametric factors $\\{\beta_{k,i}\\}$ are endowed with an
ordered spike-and-slab prior distribution. The proposed prior is critical for
coherence of the infinite-dimensional basis expansion in (1): it provides
much-needed regularization, encourages selection of a finite number of
factors, and provides posterior inference for the effective number of
nonparametric terms—including an assessment of whether any nonparametric
component is needed at all. The prior admits a parameter expansion that offers
substantial improvements in MCMC efficiency and satisfies several key
properties that broaden applicability beyond the SFFM.

The SFFM in (1) is accompanied by an observation error equation to accommodate
noisy and sparsely- or irregularly-sampled functional data. The observed data
$\bm{y}_{i}=(y_{i,1},\ldots,y_{i,m_{i}})^{\prime}$ are modeled as noisy
realizations of $Y_{i}$ on a discrete set of points
$\\{\tau_{i,j}\\}_{j=1}^{m_{i}}\subset\mathcal{T}$ for $i=1,\ldots,n$:

$y_{i,j}=Y_{i}(\tau_{i,j})+\epsilon_{i,j},\quad\epsilon_{i,j}\stackrel{{\scriptstyle
indep}}{{\sim}}N(0,\sigma_{\epsilon}^{2}),$ (2)

although non-Gaussian versions are available. We proceed using common
observation points $\tau_{i,j}=\tau_{j}$ and $m_{i}=m$ for notational
simplicity, but this restriction may be relaxed.

Although we focus on the parametric and semiparametric versions of (1), the
proposed modeling framework remains useful without any parametric template
($L=0$). In this case, (1) resembles a Karhunen-Loève decomposition, and
$\\{f_{k}\\}$ correspond to the eigenfunctions of the covariance function of
$\\{Y_{i}\\}$ (assuming the functions $Y_{i}$ have been centered). The
Karhunen-Loève decomposition provides the theoretical foundation for
functional principal components analysis (FPCA). As such, model (1)—together
with the model for $\\{f_{k}\\}$ and the ordered spike-and-slab prior for
$\\{\beta_{k,i}\\}$—constitutes a new Bayesian approach FPCA.

###  1.3 Review of related approaches

Semiparametric models for functional data are predominantly non-Bayesian. For
modeling a single function, L-splines combine a goodness-of-fit criterion with
a penalty on deviations from $\mathcal{H}_{0}$ (Ramsay and Dalzell,, 1991;
Heckman and Ramsay,, 2000). Although L-splines elegantly admit analytic
expressions for $\\{f_{k}\\}$, these derivations are highly challenging for
all but the simplest choices of $\mathcal{H}_{0}$, which inhibits widespread
practical use. In addition, L-splines offer limited direct inference on the
adequacy of the parametric templates $\mathcal{H}_{0}$, and typically require
hypothesis tests with asymptotic validity. Other non-Bayesian approaches for
semiparametric functional data analysis seek to replace nonparametric
functions with parametric alternatives. Sang et al., (2017) attempt to
simplify FPCA by using polynomials for each FPC instead of splines or Fourier
functions. In functional regression analysis, Chen et al., (2019) develop
hypothesis tests to determine whether an unknown regression function deviates
from a parametric template.

Bayesian semiparametric functional data models are less common. The
“semiparametric” model of Lee et al., (2018) refers to additive rather than
linear effects, but does not include a parametric template like
$\mathcal{H}_{0}$. Scarpa and Dunson, (2009) construct a Dirichlet process
mixture of a parametric function and a Gaussian process contamination, which
is generalized by Scarpa and Dunson, (2014) to include prior information on
the frequencies of certain functional features. These methods are designed
primarily for clustering: they identify individual curves $Y_{i}$ that deviate
substantially from the parametric model, while the remaining curves are
presumably well modeled parametrically. The SFFM is capable of modeling total
deviations from $\mathcal{H}_{0}$ for a particular $Y_{i}$, but also
captures—and corrects—partial deviations from $\mathcal{H}_{0}$ that persist
for some or many $Y_{i}$. Unlike the mixture models, the SFFM is well suited
for including additional layers of dependence, such as hierarchical (Section
6.1) or dynamic (Section 6.2) models, while maintaining efficient posterior
computing. Notably, these existing methods do not address functional
confounding.

The remainder of the paper is organized as follows. The model for the
parametric and nonparametric functions is in Section 2. The ordered spike-and-
slab prior is introduced and studied in Section 3. The MCMC algorithm is
discussed in Section 4. A simulation study is in Section 5. The model is
applied to real datasets in Section 6. We conclude in Section 7. Online
supplementary material includes proofs of all results, the full MCMC
algorithm, additional simulations, and R code for reproducibility.

##  2 Modeling the nonparametric functions

The error-free latent functions $\\{Y_{i}\\}$ belong to the space spanned by
the template parametric curves $\\{g_{\ell}\\}$ and the nonparametric curves
$\\{f_{k}\\}$. Any systemic bias resulting from the inadequacies of
$\\{g_{\ell}\\}$ must be corrected by $\\{f_{k}\\}$, which demands substantial
flexibility of the nonparametric basis $\\{f_{k}\\}$. However,
interpretability of the parametric terms $\\{g_{\ell}\\}$ and
$\\{\alpha_{\ell,i}\\}$ requires a strict distinction between the parametric
and nonparametric components.. In conjunction, these points demand both
_flexibility_ and _restraint_ from the nonparametric term.

To emphasize the importance of the distinction between $\\{g_{\ell}\\}$ and
$\\{f_{k}\\}$, consider a seemingly reasonable alternative to (1) that
replaces the nonparametric term with a Gaussian process (PFFM+gp):

$Y_{i}(\tau)=\sum_{\ell=1}^{L}\alpha_{\ell,i}g_{\ell}(\tau;\gamma)+h_{i}(\tau),\quad
h_{i}\sim\mathcal{GP}(0,\mathcal{K}_{h}).$ (3)

The PFFM+gp is a special case of the SFFM (1), and implies that each latent
curve $Y_{i}$ is a Gaussian process centered at the parametric term.
Naturally, it may be expected that each $h_{i}$ corrects for the biases of the
parametric component. However, there is a significant cost: despite being
centered at zero, the Gaussian process $h_{i}$ is not constrained to be
sufficiently distinct from the parametric term, and therefore introduces a
_functional confounding_ that biases the parametric factors
$\\{\alpha_{\ell,i}\\}$.

To illustrate this point, we fit the PFFM, PFFM+gp, and SFFM to the pinch
force data (see Section 6.1) with identical parametric models (see (11)). The
resulting posterior expectations of each $\\{\alpha_{\ell,i}\\}$ are presented
in Figure 2. Most striking, the point estimates of the parametric factors are
nearly identical between the PFFM and SFFM, while the PFFM+gp estimates are
substantially different. The discrepancy of the PFFM+gp is concerning: under
mild conditions, the PFFM will produce (approximately) unbiased estimates of
the linear coefficients $\\{\alpha_{\ell,i}\\}$ even in the presence of
dependent errors, so the point estimates are a reliable benchmark. These
results are also confirmed in the simulation study using the ground truth
values of $\\{\alpha_{\ell,i}\\}$ (Section 5). We emphasize that despite this
agreement between the PFFM and the SFFM for point estimation of
$\\{\alpha_{\ell,i}\\}$, the SFFM _does_ identify inadequacies of the PFFM and
provides key improvements in bias removal for $\\{Y_{i}\\}$ and uncertainty
quantification for both $\\{Y_{i}\\}$ and $\\{\alpha_{\ell,i}\\}$ (see
Sections 5-6).

Figure 2: Posterior expectations of $\\{\alpha_{\ell,i}\\}_{i=1}^{n}$ for
$\ell=1$ (left) and $\ell=2$ (right) for the pinch force data (see Section
6.1) comparing the PFFM estimates (x-axis) against both the SFFM and PFFM+gp
(y-axis). The point estimates of the parametric factors are nearly identical
between the PFFM and SFFM, while the PFFM+gp estimates are substantially
different. The dotted line denotes $y=x$.

The simple augmentation of a Gaussian process in (3) offers flexibility but
lacks restraint: $\\{h_{i}\\}$ absorbs variability otherwise explained by
$\\{g_{\ell}\\}$, which ultimately corrupts inference and interpretation of
the parametric term. This functional confounding is not simply an artifact of
the PFFM+gp sampling algorithm, which uses a _joint_ sampler for
$\\{\alpha_{\ell,i},h_{i}\\}$ (see the supplement), and is robust to the
choice of the covariance function $\mathcal{K}_{h}$. This issue is related to
spatial confounding (e.g., Reich et al.,, 2006) but has not been thoroughly
explored or resolved for functional data.

The proposed SFMM achieves both flexibility and distinctness by (i) modeling
each $f_{k}$ in (1) as a smooth unknown function and (ii) constraining each
$f_{k}$ to be orthogonal to $\\{g_{\ell}\\}$. The model for $f_{k}$ can be any
Bayesian curve-fitting model, such as splines, wavelets, or Gaussian
processes, usually with a prior that encourages smoothness; our specifications
are discussed subsequently. Crucially, our choice of constraints not only
ensure nonredundancy of $\\{f_{k}\\}$, but also offer key computational
simplifications that improve scalability and increase MCMC efficiency. These
results maintain regardless of the specification of $\\{g_{\ell}\\}$ or
$\\{f_{k}\\}$.

Consider model (1)-(2) evaluated at the observation points
$\\{\tau_{j}\\}_{j=1}^{m}$:

$\bm{y}_{i}=\bm{G}_{\gamma}\bm{\alpha}_{i}+\bm{F}\bm{\beta}_{i}+\bm{\epsilon}_{i},\quad\bm{\epsilon}_{i}\stackrel{{\scriptstyle
indep}}{{\sim}}N(\bm{0},\sigma_{\epsilon}^{2}\bm{I}_{m}),\quad i=1,\ldots,n$
(4)

where $\bm{G}_{\gamma}=(\bm{g}_{1;\gamma},\ldots,\bm{g}_{L;\gamma})$ is the
$m\times L$ parametric basis matrix with
$\bm{g}_{\ell;\gamma}=(g_{\ell}(\tau_{1};\gamma),\ldots,g_{\ell}(\tau_{m};\gamma))^{\prime}$,
$\bm{\alpha}_{i}=(\alpha_{1,i},\ldots,\alpha_{L,i})^{\prime}$ is the vector of
parametric factors, $\bm{F}=(\bm{f}_{1},\bm{f}_{2},\ldots)$ is the
$m\times\infty$ nonparametric basis matrix with
$\bm{f}_{k}=(f_{k}(\tau_{1}),\ldots,f_{k}(\tau_{m}))^{\prime}$, and
$\bm{\beta}_{i}=(\beta_{1,i},\beta_{2,i},\ldots,)^{\prime}$ is the vector of
nonparametric factors. By constraining each $\bm{f}_{k}$ to be orthogonal to
each $\bm{g}_{\ell;\gamma}$, the likelihood (4) factors into parametric and
nonparametric terms:

######  Lemma 1.

When $\bm{G}_{\gamma}^{\prime}\bm{f}_{k}=\bm{0}_{L}$ for all $k=1,2,\ldots$
and conditional on $\sigma_{\epsilon}^{2}$, the likelihood (4) factorizes:
$p(\bm{y}\mid\gamma,\\{\bm{\alpha}_{i}\\},\bm{F},\\{\bm{\beta}_{i}\\})=p_{0}(\bm{y}\mid\gamma,\\{\bm{\alpha}_{i}\\})p_{1}(\bm{y}\mid\bm{F},\\{\bm{\beta}_{i}\\})$,
where $p_{0}$ depends _only_ on the parametric terms $\gamma$ and
$\\{\bm{\alpha}_{i}\\}$ and $p_{1}$ depends _only_ on the nonparametric terms
$\bm{F}$ and $\\{\bm{\beta}_{i}\\}$.

Within a Bayesian model, a reasonable notion of distinctness between
parametric and nonparametric components is (conditional) independence in the
posterior. By assuming independence in the prior, Lemma 1 ensures this result:

######  Corollary 1.

Under model (4) with prior independence
$p(\\{\bm{\alpha}_{i}\\},\\{\bm{\beta}_{i}\\})=p(\\{\bm{\alpha}_{i}\\})p(\\{\bm{\beta}_{i}\\})$
and $\bm{G}_{\gamma}^{\prime}\bm{f}_{k}=\bm{0}_{L}$ for all $k=1,2,\ldots$,
the parametric and nonparametric factors are _a posteriori_ independent:
$p(\\{\bm{\alpha}_{i}\\},\\{\bm{\beta}_{i}\\}|\bm{y},\sigma_{\epsilon}^{2},\bm{F},\gamma)=p(\\{\bm{\alpha}_{i}\\}|\bm{y},\sigma_{\epsilon}^{2},\bm{F},\gamma)p(\\{\bm{\beta}_{i}\\}|\bm{y},\sigma_{\epsilon}^{2},\bm{F},\gamma)$.

Corollary 1 is crucial not only for preserving _distinctness_ between the
parametric and nonparametric terms, but also for designing an _efficient_ MCMC
sampler. The specific constraint
$\bm{G}_{\gamma}^{\prime}\bm{f}_{k}=\bm{0}_{L}$ ensures that the sampling
steps for all parametric factors $\\{\bm{\alpha}_{i}\\}$ do _not_ depend on
the nonparametric parameters in any way: the sampling steps for
$\\{\bm{\alpha}_{i}\\}$ are identical to the fully parametric case, while the
sampling steps for $\\{\bm{\beta}_{i}\\}$ proceed exactly as in a fully
nonparametric setting (see Section 3). This decoupling of parametric and
nonparametric sampling steps is particularly useful for more complex models.
At the same time, Corollary 1 implies that _separate_ sampling steps for
$\\{\bm{\alpha}_{i}\\}$ and $\\{\bm{\beta}_{i}\\}$ are equivalent to a _joint_
sampler for all $\\{\bm{\alpha}_{i},\bm{\beta}_{i}\\}$. Joint sampling steps
are typically preferred for greater MCMC efficiency, but often increase the
computational burden. For instance, a joint sampler for
$\\{\bm{\alpha}_{i},h_{i}\\}$ under the PFFM+gp requires either a large block
sampler or marginalization over either $\bm{\alpha}_{i}$ or $h_{i}$, both of
which become more challenging when the model for $\bm{\alpha}_{i}$ is complex.
The SFFM—with the proposed orthogonality constraint—completely avoids this
tradeoff and guarantees joint sampling steps even for complex models. Such a
reduction is not available in general, for example $L^{2}$ orthogonality
$\int_{\mathcal{T}}g_{\ell}(\tau;\gamma)f_{k}(\tau)\,d\tau=0$ for
$\ell=1,\ldots,L$.

To enforce the proposed constraint, we _condition_ on orthogonality
$\bm{G}_{\gamma}^{\prime}\bm{f}_{k}=\bm{0}_{L}$ in the model for $\bm{f}_{k}$.
Conditioning is a natural Bayesian mechanism for incorporating such
information into estimation and inference. We model the unknown functions
$f_{k}(\tau)=\bm{b}^{\prime}(\tau)\bm{\psi}_{k}$ using known basis functions
$\bm{b}$ and unknown coefficients $\bm{\psi}_{k}$, so a prior on
$\bm{\psi}_{k}$ implies a prior on the function $f_{k}$, for example to
encourage smoothness. Our default specification uses low-rank thin plate
splines $\bm{b}$, which are flexible, computationally efficient, and well-
defined on $\mathcal{T}\subset\mathbb{R}^{d}$ for $d\in\mathbb{Z}^{+}$, along
with the smoothness prior
$[\bm{\psi}_{k}\mid\lambda_{f_{k}}]\stackrel{{\scriptstyle
indep}}{{\sim}}N(\bm{0},\lambda_{f_{k}}^{-1}\bm{\Omega}^{-})$, where
$\bm{\Omega}$ is the (known) matrix of integrated squared second derivatives,
$[\bm{\Omega}]_{j,j^{\prime}}=\int\ddot{b}_{j}(\tau)\ddot{b}_{j^{\prime}}(\tau)d\tau$,
which implies that the prior distribution represents the classical roughness
penalty $-2\log p(\bm{\psi}_{k}\mid\lambda_{f_{k}})\stackrel{{\scriptstyle
c}}{{=}}\lambda_{f_{k}}\bm{\psi}_{k}^{\prime}\bm{\Omega}\bm{\psi}_{k}=\lambda_{f_{k}}\int_{\mathcal{T}}\\{\ddot{f}_{k}(\tau)\\}^{2}d\tau$
and the smoothing parameter $\lambda_{f_{k}}^{-1/2}\stackrel{{\scriptstyle
indep}}{{\sim}}\mbox{Uniform}(0,10^{4})$ is learned from the data.

Absent the orthogonality constraint, the full conditional distribution of each
vector of coefficients is $[\bm{\psi}_{k}\mid-]\sim
N\left(\bm{Q}_{\psi_{k}}^{-1}\bm{\ell}_{\psi_{k}},\bm{Q}_{\psi_{k}}^{-1}\right)$,
where
$\bm{Q}_{\psi_{k}}=\sigma_{\epsilon}^{-2}(\bm{B}^{\prime}\bm{B})\sum_{i=1}^{n}\beta_{k,i}^{2}+\lambda_{f_{k}}\bm{\Omega}$,
$\bm{\ell}_{\psi_{k}}=\sigma_{\epsilon}^{-2}\bm{B}^{\prime}\sum_{i=1}^{n}\big{\\{}\beta_{k,i}\big{(}\bm{y}_{i}-\bm{G}_{\gamma}\bm{\alpha}_{i}-\sum_{\ell\neq
k}\bm{f}_{\ell}\beta_{\ell,i}\big{)}\big{\\}}$, and
$\bm{B}=(\bm{b}^{\prime}(\tau_{1}),\ldots,\bm{b}^{\prime}(\tau_{m}))^{\prime}$.
Samples from the full conditional distribution of $\bm{\psi}_{k}$
_conditional_ on $\bm{G}_{\gamma}^{\prime}\bm{f}_{k}=\bm{0}_{L}$ can be
obtained by (i) sampling from the unconstrained distribution
$[\bm{\psi}_{k}\mid-]$ and (ii) applying a simple shift:

######  Lemma 2.

Suppose $\bm{G}_{\gamma}$ is $L\times m$ with rank $L$. Denote
$\bm{\psi}_{k}^{0}\sim
N\left(\bm{Q}_{\psi_{k}}^{-1}\bm{\ell}_{\psi_{k}},\bm{Q}_{\psi_{k}}^{-1}\right)$
with $\bm{f}_{k}^{0}=\bm{B}\bm{\psi}_{k}^{0}$. The shifted term
$\bm{f}_{k}=\bm{B}\bm{\psi}_{k}$ with
$\bm{\psi}_{k}=\bm{\psi}_{k}^{0}-\bm{Q}_{\psi_{k}}^{-1}\bm{B}^{\prime}\bm{G}_{\gamma}\left(\bm{G}_{\gamma}^{\prime}\bm{B}\bm{Q}_{\psi_{k}}^{-1}\bm{B}^{\prime}\bm{G}_{\gamma}\right)^{-1}\bm{G}_{\gamma}^{\prime}\bm{B}\bm{\psi}_{k}^{0}$
satisfies $\bm{f}_{k}\stackrel{{\scriptstyle
d}}{{=}}[\bm{f}_{k}^{0}\mid\bm{G}_{\gamma}^{\prime}\bm{f}_{k}^{0}=\bm{0}]$ and
$\mathbb{P}(\bm{G}_{\gamma}^{\prime}\bm{f}_{k}=\bm{0})=1$.

Lemma 2 provides a simple and effective mechanism to adapt any standard
(basis) regression sampling step for $[\bm{\psi}_{k}\mid-]$ to accommodate
this crucial orthogonality constraint, including in the case when the
nonlinear term $\gamma$ is unknown. Note that it is possible to adapt Lemma 2
for the PFFM+gp; however, the resulting model still does not provide direct
inference on the necessity of the nonparametric factors (i.e., rank selection)
or describe the systemic biases of $\mathcal{H}_{0}$, which are captured by
$\\{f_{k}\\}$ in the SFFM. With or without additional constraints, the PFFM+gp
sampling algorithm lacks the same scalability in $m$ (see Section 6).

For further computational simplifications along with parameter
identifiability, we additionally constrain $\bm{G}_{\gamma}$ and $\bm{F}$ to
be _orthonormal_ , which implies that the joint basis
$\\{g_{1},\ldots,g_{L},f_{1},f_{2},\ldots\\}$ in (1) is orthonormal. First,
$\bm{G}_{\gamma}$ is orthogonalized using a QR decomposition. For any value of
$\gamma$, let $\bm{G}_{\gamma}^{0}=\bm{Q}_{\gamma}\bm{R}_{\gamma}$ be the QR
decomposition of the initial basis matrix
$\bm{G}_{\gamma}^{0}=(\bm{g}_{1;\gamma},\ldots,\bm{g}_{L;\gamma})$. By setting
$\bm{G}_{\gamma}=\bm{Q}_{\gamma}$, we ensure that
$\bm{G}_{\gamma}^{\prime}\bm{G}_{\gamma}=\bm{I}_{L}$ and the columns of
$\bm{G}_{\gamma}$ span the same space as the columns of $\bm{G}_{\gamma}^{0}$.
When $\gamma$ is unknown and endowed with a prior distribution, the QR
decomposition is incorporated into the likelihood evaluations of (4) for
posterior sampling of $\gamma$. The orthonormality of $\bm{G}_{\gamma}$
provides critical simplifications for efficient posterior inference:

######  Corollary 2.

Under model (4) and subject to $\bm{G}_{\gamma}^{\prime}\bm{f}_{k}=\bm{0}_{L}$
for all $k$ and $\bm{G}_{\gamma}^{\prime}\bm{G}_{\gamma}=\bm{I}_{L}$, the
likelihood for the parametric factors is proportional to
$p(\bm{y}\mid\gamma,\\{\bm{\alpha}_{i}\\},\bm{F},\\{\bm{\beta}_{i}\\})\propto
p_{0,G}(\bm{y}\mid\gamma,\\{\bm{\alpha}_{i}\\})$ where $p_{0,G}$ is the
likelihood defined by
$[\bm{g}_{\ell;\gamma}^{\prime}\bm{y}_{i}\mid-]\stackrel{{\scriptstyle
indep}}{{\sim}}N(\alpha_{\ell,i},\sigma_{\epsilon}^{2})$ for
$\ell=1,\ldots,L$.

Corollary 2 implies that posterior inference for the parametric factors only
requires a simple independent Gaussian likelihood with pseudo-data
$\\{\bm{g}_{\ell;\gamma}^{\prime}\bm{y}_{i}\\}$. Hence, the posterior sampling
steps for $\\{\alpha_{\ell,i}\\}$ are often convenient: for example, under the
prior $\alpha_{\ell,i}\stackrel{{\scriptstyle
indep}}{{\sim}}N(0,\sigma_{\alpha_{\ell}}^{2})$, the joint full conditional
distribution of $\\{\alpha_{\ell_{i}}\\}$ is decomposable into independent
univariate Gaussian distributions with
$\mathbb{E}(\alpha_{\ell,i}\mid-)=\sigma_{\alpha_{\ell}}^{2}/(\sigma_{\alpha_{\ell}}^{2}+\sigma_{\epsilon}^{2})\bm{g}_{\ell;\gamma}^{\prime}\bm{y}_{i}$
and
$\mbox{Var}(\alpha_{\ell,i}\mid-)=\sigma_{\epsilon}^{2}\sigma_{\alpha_{\ell}}^{2}/(\sigma_{\alpha_{\ell}}^{2}+\sigma_{\epsilon}^{2})$,
which may be sampled efficiently—and is identical for both the PFFM and the
SFFM. Similar simplifications apply for more complex (e.g., dynamic) models
for the parametric factors, and are invariant to the parametrization of the
remaining model. We emphasize that the constraint on $\bm{G}_{\gamma}$ is
simply a computationally favorable reparametrization of the parametric term
and does not hinder the interpretation: since
$\bm{G}_{\gamma}^{0}\bm{\alpha}_{i}^{0}=\bm{G}_{\gamma}\bm{\alpha}_{i}$ for
$\bm{\alpha}_{i}=\bm{R}_{\gamma}\bm{\alpha}_{i}^{0}$ under the QR
decomposition, we recover the parametric factors on the original scale by
setting $\bm{\alpha}_{i}^{0}=\bm{R}_{\gamma}^{-1}\bm{\alpha}_{i}$, which can
be computed draw by draw within the MCMC sampler.

For each nonparametric term $\bm{f}_{k}$, we incorporate the _orthogonality_
constraints via conditioning on $\bm{f}_{k^{\prime}}^{\prime}\bm{f}_{k}=0$ for
all $k^{\prime}\neq k$ and the _unit-norm_ constraint by rescaling
$\bm{f}_{k}$ appropriately. These steps follow Kowal, (2020) and are described
briefly. The linear orthogonality constraint is enforced using a
straightforward modification of Lemma 2, which augments the parametric
orthogonality of $\bm{G}_{\gamma}$ with the nonparametric terms
$\bm{f}_{k^{\prime}}$ for $k^{\prime}\neq k$. The unit-norm constraint is
enforced by suitably rescaling $\bm{\psi}_{k}$ after sampling; we equivalently
rescale $\beta_{k,i}$ to preserve the the product $\bm{f}_{k}\beta_{k,i}$ in
the likelihood (4). This operation does not change the shape of the curve
$f_{k}$ nor the likelihood (4).

Naturally, Corollary 2 may be adapted for the nonparametric factors:

######  Corollary 3.

Under model (4) and subject to $\bm{G}_{\gamma}^{\prime}\bm{f}_{k}=\bm{0}_{L}$
for all $k$ and $\bm{F}^{\prime}\bm{F}=\bm{I}_{\infty}$, the likelihood for
the nonparametric factors is proportional to
$p(\bm{y}\mid\gamma,\\{\bm{\alpha}_{i}\\},\bm{F},\\{\bm{\beta}_{i}\\})\propto
p_{1,F}(\bm{y}\mid\bm{F},\\{\bm{\beta}_{i}\\})$ where $p_{1,F}$ is the
likelihood defined by
$[\bm{f}_{k}^{\prime}\bm{y}_{i}\mid-]\stackrel{{\scriptstyle
indep}}{{\sim}}N(\beta_{k,i},\sigma_{\epsilon}^{2})$ for all $k$.

As with the parametric factors, Corollary 3 implies that posterior inference
for the nonparametric factors only requires a simple independent Gaussian
likelihood with pseudo-data $\\{\bm{f}_{k}^{\prime}\bm{y}_{i}\\}$. As a
result, the SFFM is capable of incorporating complex prior specifications for
$\\{\beta_{k,i}\\}$—which are essential for inference in the SFFM with the
_infinite_ expansion in (1).

##  3 Modeling the nonparametric factors with rank selection

Regularization of the SFFM is critical: the infinite-dimensional nonparametric
term in (1) is clearly overparametrized for finite $m$. The nonparametric
components increase model complexity and should be removed whenever the added
complexity is not supported by the data. This _rank selection_ is complicated
by the presence of the parametric template with unknown
$\\{\alpha_{\ell,i}\\}$ and $\gamma$. Information criteria offer metrics for
model comparisons, but require separate model fits for each rank. Ordered
shrinkage priors (e.g., Bhattacharya and Dunson,, 2011) can reduce sensitivity
to the choice of the rank, but often induce undesirable overshrinkage.
Instead, we design an _ordered spike-and-slab_ prior that applies joint
shrinkage and selection to the nonparametric _factors_
$\\{\beta_{k,i}\\}_{i=1}^{n}$ and therefore removes unnecessary nonparametric
terms from the SFFM. Because of the proposed orthogonality constraints on
$\\{g_{\ell}\\}$ and $\\{f_{k}\\}$ and the accompanying simplifications from
Corollaries 2 and 3, we are able to introduce a sophisticated prior for
$\\{\beta_{k,i}\\}$ within the broader SFFM framework of (1)—and do so with
minimal impact on computational cost and MCMC efficiency relative to the PFFM
(see Section 6).

###  3.1 Ordered spike-and-slab priors

Let $\theta_{k}$ denote the parameter that determines the inclusion of the
$k$th nonparametric factor $\\{\beta_{k,i}\\}_{i=1}^{n}$; the connection is
made explicit below. We apply a spike-and-slab prior for
$\\{\theta_{k}\\}_{k=1}^{\infty}$ that imposes a particular _ordering_
structure on the spike probabilities $\\{\pi_{k}\\}_{k=1}^{\infty}$:

$\displaystyle[\theta_{k}\mid\pi_{k}]$ $\displaystyle\sim P_{k},\quad
P_{k}=(1-\pi_{k})P_{slab}+\pi_{k}P_{spike}$ (5)
$\displaystyle\pi_{k}=\sum_{h=1}^{k}\omega_{h},\quad\omega_{h}=\nu_{h}\prod_{\ell=1}^{h-1}(1-\nu_{\ell}),\quad\nu_{\ell}\stackrel{{\scriptstyle
iid}}{{\sim}}\mbox{Beta}(\iota,\iota\kappa),$ (6)

where $P_{slab}$ is distribution of the slab (active) component and
$P_{spike}$ is the distribution of the spike (inactive) component. The
cumulative summation for $\pi_{k}$ ensures an increasing sequence of (spike)
probabilities, $\pi_{k}<\pi_{k+1}$, that converges:
$\lim_{k\rightarrow\infty}\pi_{k}=1$. The hyperparameters are determined by
the distributions $P_{slab}$ and $P_{spike}$ and the scalars $\iota,\kappa>0$.
The prior in (5)-(6) generalizes the _cumulative shrinkage process_ (CUSP)
introduced by Legramanti et al., (2020), who specified
$P_{spike}=\delta_{\theta_{\infty}}$ to be a point mass at $\theta_{\infty}$
with $\iota=1$ and fixed $\kappa$ (using $\kappa=5$ in their examples). The
CUSP in (5)-(6) will be further adapted for the SFFM model (1) in Section 3.3.

The ordering of the spike probabilities $\\{\pi_{k}\\}_{k=1}^{\infty}$ in (6)
implies an ordering for the prior distribution of the parameters
$\\{\theta_{k}\\}_{k=1}^{\infty}$ in (5):

######  Proposition 1.

For $\varepsilon>0$ and fixed $\theta_{0}$, let
$\mathbb{B}_{\varepsilon}(\theta_{0})=\\{\theta_{k}:|\theta_{k}-\theta_{0}|<\varepsilon\\}$.
Prior (5)-(6) implies that
$\mathbb{P}(|\theta_{k}-\theta_{0}|\leq\varepsilon)<\mathbb{P}(|\theta_{k+1}-\theta_{0}|\leq\varepsilon)$
whenever
$P_{slab}\\{\mathbb{B}_{\varepsilon}(\theta_{0})\\}<P_{spike}\\{\mathbb{B}_{\varepsilon}(\theta_{0})\\}$.

Intuitively, for any $\theta_{0}$ that is favored under the spike distribution
$P_{spike}$ relative to the slab distribution $P_{slab}$, the CUSP prior
(5)-(6) places greater mass around $\theta_{0}$ as $k$ increases. For spike-
and-slab priors, we are most interested in $\theta_{0}=0$. The special case of
Proposition 1 with $P_{spike}=\delta_{\theta_{\infty}}$, $\iota=1$, and
$\theta_{0}=0$ is proved by Legramanti et al., (2020).

To specify $P_{slab}$ and $P_{spike}$, we apply a normal mixture of inverse-
gamma (NMIG) prior:

$[\eta_{k}|\theta_{k},\sigma_{k}^{2}]\stackrel{{\scriptstyle
indep}}{{\sim}}N(0,\theta_{k}\sigma_{k}^{2}),\quad[\theta_{k}|\pi_{k}]\sim(1-\pi_{k})\delta_{1}+\pi_{k}\delta_{v_{0}},\quad[\sigma_{k}^{-2}]\stackrel{{\scriptstyle
iid}}{{\sim}}\mbox{Gamma}(a_{1},a_{2})$ (7)

where the scalars $\\{\eta_{k}\\}$ are linked to the factors
$\\{\beta_{k,i}\\}$ in Section 3.3 and $v_{0}$, $a_{1}$ and $a_{2}$ are
hyperparameters. The NMIG prior incorporates variable selection by assigning
the variance scale parameter $\theta_{k}$ to the slab component,
$\theta_{k}=1$, or the spike component, $\theta_{k}=v_{0}$. By design, the
NMIG prior produces a continuous distribution for the conditional variance of
$\eta_{k}$, which is preferable for variable selection and risk properties
(Ishwaran and Rao,, 2005). The NMIG prior is often less sensitive to
hyperparameter choices compared to other spike-and-slab priors (Ishwaran and
Rao,, 2005; Scheipl et al.,, 2012).

By coupling the NMIG prior for $\\{\eta_{k}\\}$ in (7) with the CUSP prior for
$\\{\pi_{k}\\}$ in (6) and marginalizing over
$\\{\theta_{k},\sigma_{k}^{2},\pi_{k}\\}$, we obtain the following ordering
result for the marginal prior on $\\{\eta_{k}\\}_{k=1}^{\infty}$:

######  Corollary 4.

For $\varepsilon>0$,
$\mathbb{P}(|\eta_{k}|\leq\varepsilon)<\mathbb{P}(|\eta_{k+1}|\leq\varepsilon)$
whenever $v_{0}<1$.

Corollary 4 follows from Proposition 1 by observing that
$[\eta_{k}|\pi_{k}]\sim(1-\pi_{k})t_{2a_{1}}(0,\sqrt{a_{2}/a_{1}})+\pi_{k}t_{2a_{1}}(0,\sqrt{v_{0}a_{2}/a_{1}})$,
where $t_{d}(m,s)$ denotes a $t$-distribution with mean $m$, standard
deviation $s$, and degrees of freedom $d$. This mixture is a special case of
(5), and the densities of the $t$-distributions place greater mass near zero
when the scale parameter is smaller. Hence, the condition $v_{0}<1$ ensures
that the spike distribution is indeed more concentrated around zero.

For interpretability and efficient MCMC sampling, there is a convenient data
augmentation of the CUSP prior (5). Let $z_{k}\in\\{1,\ldots,\infty\\}$ denote
a categorical variable with $\mathbb{P}(z_{k}=h|\omega_{h})=\omega_{h}$. The
specification $[\theta_{k}|z_{k}]\sim(1-\mathbb{I}\\{z_{k}\leq
k\\})P_{slab}+\mathbb{I}\\{z_{k}\leq k\\}P_{spike}$ induces (5) via
marginalization over $z_{k}$. The number of active (slab) terms is therefore
$K^{*}=\sum_{k=1}^{\infty}\mathbb{I}\\{z_{k}>k\\},$ with posterior inference
available through the proposed MCMC sampling algorithm (see Algorithm 1). For
model (1), $K^{*}$ is the effective number of nonparametric terms, and
therefore is an important inferential target to assess the adequacy of the
parametric model.

The prior expected number of slab terms is $\mathbb{E}(K^{*})=\kappa$, which
appears in (6) and allows for the inclusion of prior information regarding the
number of factors. The choice of $\kappa$ can influence the posterior
distribution for $K^{*}$, and thus subsequent inference on the number of
nonparametric factors in (1). To mitigate this effect, we propose the
hyperprior $\kappa\sim\mbox{Gamma}(a_{\kappa},b_{\kappa})$ for
$a_{\kappa},b_{\kappa}>0$, which is conditionally conjugate to (6). The
hyperparameters may be selected to provide weak prior information regarding
the number of factors: in practice, we set $a_{\kappa}=2$ and $b_{\kappa}=1$
so that $\mathbb{E}(\kappa)=\mbox{Var}(\kappa)=2$. Similarly, we select the
default value $\iota=1$ for simplicity.

CUSPs are intrinsically linked to Indian buffet processes (IBPs; Griffiths and
Ghahramani,, 2011), which are often used for sparse (non-functional) factor
models with rank selection (Ohn and Kim,, 2021). IBPs supply a prior over
binary matrices $\\{b_{j,k}\\}$ with $m$ rows and infinitely many columns, and
are usually applied to $\bm{F}$ in non-functional versions of (4). IBPs are
obtained by establishing the conditional distributions
$[b_{j,k}\mid\mu_{k}]\stackrel{{\scriptstyle
indep}}{{\sim}}\mbox{Bernoulli}(\mu_{k})$ and
$[\mu_{k}]\stackrel{{\scriptstyle
iid}}{{\sim}}\mbox{Beta}(\iota\kappa/K,\iota)$, usually with $\iota=1$, and
then integrating over $\\{\mu_{k}\\}$ with $K\rightarrow\infty$. IBPs admit a
stick-breaking construction (Teh et al.,, 2007), which can be mapped to the
CUSP probability sequence $\\{\pi_{k}\\}_{k=1}^{\infty}$:

######  Proposition 2.

The CUSP (6) satisfies $(1-\pi_{k})=\mu_{(k)}$ where
$\mu_{(1)}>\cdots>\mu_{(K)}$ are the ordered (slab) probabilities from the Teh
et al., (2007) stick-breaking construction of the IBP, i.e.,
$\mu_{(k)}=\prod_{\ell=1}^{k}\nu_{\ell}^{\prime}$ with
$\nu_{\ell}^{\prime}\stackrel{{\scriptstyle
iid}}{{\sim}}\mbox{Beta}(\iota\kappa,\iota).$

Conventionally, IBPs apply the multiplicative beta process to the slab
probability $(1-\pi_{k})$, which explains the complement in Proposition 2.

Despite these fundamental stick-breaking connections, CUSPs and IBPs remain
distinct due to (5): CUSPs define a spike-and-slab prior for a sequence of
parameters $\\{\theta_{k}\\}_{k=1}^{\infty}$, while IBPs define a prior over
binary matrices. In the context of (1), elementwise sparsity of $\bm{F}$ is
unwarranted: although we are interested in rank selection, we prefer priors
that encourage smoothness rather than sparsity in $f_{k}(\tau)$ over
$\tau\in\mathcal{T}$.

###  3.2 Consistent rank selection

To motivate the use of this prior within the SFFM, we investigate the
asymptotic behavior of the generalized CUSP as a standalone prior for rank
selection. Specifically, we consider the modified setting with

$y_{i}=\eta_{0i}+\epsilon_{i},\quad{\epsilon_{i}}\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,1),$ (8)

where $\eta_{01},\dots,\eta_{0n}$ are “true” mean parameters. This setting is
applicable to the SFFM via Corollary 3 and is related to the problem of
estimating a high-dimensional mean vector from a single multivariate
observation (Castillo and van der Vaart,, 2012; Rockova,, 2018), but with
ordered sparsity $\eta_{0k}\neq 0$ for $k\leq K_{0n}$ and $\eta_{0k}=0$ for
$k>K_{0n}$. As such, $K_{0n}$ is the true number of nonzero means and
corresponds to the true rank in model (1). As customary in posterior
asymptotics, we allow $K_{0n}$ to grow with $n$.

We study the posterior of $\eta_{i}$ under the CUSP (5)-(6) and NMIG (7)
prior. A key term is the “remainder” $R_{n}=\sum_{k\geq K_{0n}}\omega_{k}$,
which represents an upper bound for the probability mass assigned to the slab
for the parameters that are null. Under the CUSP, we show that $R_{n}$ is far
from zero with vanishing probability _a priori_ :

######  Lemma 3.

Let $\varepsilon_{n}\to 0$ with
$\varepsilon_{n}^{1/K_{0n}}>\kappa/(\kappa+1)$. For the CUSP prior (5)-(6) and
a positive constant $C>1$, the remainder term $R_{n}=\sum_{k\geq
K_{0n}}\omega_{k}$ satisfies
$\mathbb{P}(R_{n}>\varepsilon_{n})\leq\exp(-CK_{0n}).$

According to Lemma 6, the prior probability that $R_{n}$ exceeds a small
threshold is exponentially small. Such exponential decay in the prior is
essential to obtain optimal posterior behavior in sparse settings (Castillo
and van der Vaart,, 2012). This result does not require the NMIG prior (7).

To connect the prior behavior in Lemma 6 to the posterior, first consider the
following mild conditions on the prior parameters: (C1)
$a_{1}/(a_{2}v_{0})>n^{2}/K_{0n}$; (C2) $AK_{0n}/n<1/2$ with $A>1/2$; (C3)
$a_{1}>1/2$, and $a_{1}\leq a_{2}$. Condition (C1) implies that the precision
of the spike component is increasing with $n$; condition (C2) requires that
$K_{0n}$ is strictly less than the total number of elements (or the truncation
limit); and condition (C3) implies that the values of the $\sigma_{k}^{2}$
parameters in (7) are not too large. Now let $\mathbb{P}_{0}$ and
$\mathbb{E}_{0}$ denote the probability and the expectation, respectively,
under the true distribution of the data (8). Under (C1)-(C3), the posterior
probability that the remainder exceeds a shrinking $\varepsilon_{n}$, for
$n\to\infty$, goes to zero in expectation $\mathbb{E}_{0}$. This is formalized
in the following result:

######  Theorem 1.

Let $\varepsilon_{n}\to 0$ with $\varepsilon_{n}^{1/K_{0n}}>\kappa/(\kappa+1)$
and assume (C1)-(C3) and $C>2Ae$. For the CUSP (5)-(6) and NMIG (7) priors,
the posterior distribution satisfies

$\lim_{n\to\infty}\mathbb{E}_{0}\left\\{\mathbb{P}\left(R_{n}>\varepsilon_{n}\mid
y_{1},\ldots,y_{n}\right)\right\\}=0.$

In the context of the SFFM, (8) can be obtained by fixing $\bm{F}$ and setting
$\beta_{k,i}=\eta_{k,i}$. For simplicity, set the observation error variance
in (1) to be $\sigma_{\epsilon}^{2}=n$, which grows with the number of curves.
Applying Corollary 3, the likelihood for $\\{\eta_{k}\\}$ in (8) is equivalent
to the likelihood implied by $y_{k}=\bar{\eta}_{k}+\epsilon_{k}$ with
$\epsilon_{k}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,1)$, where
$y_{k}=n^{-1}\sum_{i=1}^{n}\bm{f}_{k}^{\prime}\bm{y}_{i}$ and
$\bar{\eta}_{k}=n^{-1}\sum_{i=1}^{n}\eta_{k,i}$. The above results then apply
for the NMIG prior on $\\{\bar{\eta}_{k}\\}$ and by interchanging the notation
$k$ for $n$.

###  3.3 Parameter-expanded spike-and-slab models

The generalized CUSP provides an appealing mechanism for introducing ordered
shrinkage and selection into model (1). A subtle yet important feature of this
framework is that spike-and-slab selection is applied jointly to all $n$
factors $\\{\beta_{k,i}\\}_{i=1}^{n}$. However, parametrization of a
_collection_ of variables via a spike-and-slab prior on a _scalar_ $\eta_{k}$
is known to introduce significant challenges for MCMC sampling. In a related
setting, Scheipl et al., (2012) showed that as the size of this collection
increases (here, $n$), it becomes more difficult for the MCMC algorithm to
switch between the slab and spike components. The consequences are nontrivial:
without the ability to traverse between the slab and spike components, the
uncertainty quantification for each factor’s activeness is invalidated, and
inference for the effective number of nonparametric factors $K^{*}$ is not
reliable.

To illustrate the implications, we implemented the CUSP prior and sampling
algorithm from Legramanti et al., (2020), which does not address these issues.
The analogous prior in Legramanti et al., (2020) is
$[\beta_{k,i}\mid\theta_{k}]\stackrel{{\scriptstyle
indep}}{{\sim}}N(0,\theta_{k})$ and a CUSP prior for
$\\{\theta_{k}\\}_{k=1}^{\infty}$ with $P_{spike}=\delta_{\theta_{\infty}}$
and $P_{slab}=\mbox{Inv-Gamma}(a_{\theta},b_{\theta})$. We use the
hyperparameters from Legramanti et al., (2020): $\theta_{\infty}=0.05$,
$a_{\theta}=b_{\theta}=2$, and $\kappa=5$, but the subsequent results persist
for other choices. For this example, we use the simulation design from Section
5 with $K_{true}=6$ nonparametric factors. Figure 3 (left) displays the MCMC
samples of $K^{*}$ across five different MCMC runs with random
initializations. The Legramanti et al., (2020) prior and sampler produces MCMC
draws that are essentially stuck at a single value: $K^{*}$ changes values in
only 0.004% of the MCMC iterations. Clearly, this result cannot be used to
provide reliable estimation or uncertainty quantification for the effective
number of nonparametric components $K^{*}$.

Figure 3: MCMC samples of $K^{*}$ across 5 chains for Legramanti et al.,
(2020) (left) compared to the proposed prior and algorithm (right). The
horizontal black line denotes the true value $K_{true}=6$. The proposed
approach provides substantial improvements in mixing and produces more
accurate estimates.

Our solution is to introduce a redundant parameter expansion for the ordered
spike-and-slab prior. We induce a prior on each factor $\beta_{k,i}$ via the
following parameter expansion:

$\beta_{k,i}=\eta_{k}\xi_{k,i},\quad[\xi_{k,i}\mid
m_{\xi_{k,i}}]\stackrel{{\scriptstyle indep}}{{\sim}}N(m_{\xi_{k,i}},1),\quad
m_{\xi_{k,i}}\stackrel{{\scriptstyle
indep}}{{\sim}}\frac{1}{2}\delta_{1}+\frac{1}{2}\delta_{-1}$ (9)

where $\eta_{k}$ applies the ordered spike-and-slab shrinkage for factor $k$
via the NMIG prior (7) and $\\{\xi_{k,i}\\}_{i=1}^{n}$ disperses this
shrinkage throughout $\\{\beta_{k,i}\\}_{i=1}^{n}$. Here, $\xi_{k,i}$ is
nonidentifiable and exists primarily for the purpose of improving MCMC
performance. The Gaussian parameter expansion of $\xi_{k,i}$ is centered
around $1$ or $-1$ with variance $1$. As a result, the parameter $\eta_{k}$
retains the interpretability as a selection parameter and assigns the
magnitude to $\\{\beta_{k,i}\\}_{i=1}^{n}$.

The impact of the parameter expansion is striking: Figure 3 (right) shows the
MCMC samples of $K^{*}$ from the proposed method (again across five different
MCMC runs with random initializations). The proposed approach produces draws
that are centered at the true value with some variability: $K^{*}$ changes
values in 32% of the MCMC iterations, which is a significant improvement over
Legramanti et al., (2020). Clearly, the MCMC-based uncertainty quantification
from the proposed approach is more plausible—which is essential for inference
on effective number of nonparametric terms in the SFFM (1).

The parameter-expanded prior (9) also maintains the required ordering on the
factors:

######  Proposition 3.

For $\varepsilon>0$,
$\mathbb{P}(|\beta_{k,i}|\leq\varepsilon)<\mathbb{P}(|\beta_{k+1,i}|\leq\varepsilon)$
whenever $v_{0}<1$.

Proposition 3 requires only that $\\{\xi_{k,i}\\}$ are marginally independent
and identically distributed. Therefore, other parameter expansion schemes may
be substituted in (9) without disrupting the probabilistic ordering of
$|\beta_{k,i}|$. By a similar argument, the ordering also is preserved across
observations:
$\mathbb{P}(|\beta_{k,i}|\leq\varepsilon)<\mathbb{P}(|\beta_{k+1,i^{\prime}}|\leq\varepsilon)$
for $\varepsilon>0$ and $v_{0}<1$.

##  4 MCMC for posterior inference

We design an efficient MCMC algorithm for posterior inference. The proposed
Gibbs sampler is both computationally scalable and MCMC efficient (see Section
6) due to two key features: the orthogonality constraints on $\\{g_{\ell}\\}$
and $\\{f_{k}\\}$, which produce important model simplifications, and the
parameter-expanded sampler for the ordered spike-and-slab prior, which
substantially improves MCMC mixing (see Figure 3). We implement the following
sampling steps, with details provided in the supplement: (i) imputation of
$y_{i}(\tau^{*})$ for any unobserved $\tau^{*}$; (ii) sampling the nonlinear
parameter $\gamma$ if unknown; (iii) sampling the constrained nonparametric
factors $\\{f_{k}\\}$ (see Section 2) along with the smoothing parameters
$\\{\lambda_{f_{k}}\\}$; (iv) sampling the ordered spike-and-slab parameters
(Algorithm 1); (v) sampling the parametric factors $\\{\alpha_{\ell,i}\\}$
using the simplifications from Corollary 2; and (vi) sampling any variance
components, such as $\sigma_{\epsilon}^{2}$.

The implementation incorporates a truncation of the infinite summation in (1),
which provides simpler and faster computations. The following result shows
that finite approximations are accurate for sufficiently large truncation $K$:

######  Proposition 4.

Let $\theta^{(K)}=\\{\theta_{k}\\}_{k=1}^{K}$ denote the sequence
$\\{\theta_{k}\\}_{k=1}^{\infty}$ truncated at $K$. Under the CUSP (5)-(6) and
NMIG (7) priors with $0<v_{0}<\varepsilon<1$, we have
$\mathbb{P}\left\\{d_{\infty}(\theta,\theta^{(K)})>\varepsilon\right\\}\leq\kappa\\{\kappa/(1+\kappa)\\}^{K}$.

The approximation error induced by truncating
$\\{\theta_{k}\\}_{k=1}^{\infty}$ to $K$ terms decreases rapidly in $K$, which
suggests that the proposed infinite-dimensional ordered spike-and-slab prior
is accurately approximated by a conservative truncation.

We emphasize that the sampling steps for the ordered spike-and-slab
parameters, detailed in Algorithm 1, consist of simple, fast, and closed form
updates. Because of the orthogonality constraints and Corollary 3, these
distributions depend on the data only through
$\\{\bm{f}_{k}^{\prime}\bm{y}_{i}\\}$. In addition, we retain posterior
samples of $K^{*}=\sum_{k=1}^{K}\mathbb{I}\\{z_{k}>k\\}$, which is the
effective number of nonparametric terms. If the posterior distribution of
$K^{*}$ places mass on values near $K$, then the truncation level $K$ should
be increased. Adaptations of this sampling algorithm for use of the ordered
spike-and-slab prior with other models, such as (non-functional) factor
models, are available by replacing $y_{k,i}^{F}$ with the appropriate term.

Let $y_{k,i}^{F}=\bm{f}_{k}^{\prime}\bm{y}_{i}$ for $i=1,\ldots,n$ and
$k=1,\ldots,K$:  1. Sample $[m_{\xi_{k,i}}|-]$ from $\\{-1,1\\}$ with
$\mathbb{P}(m_{\xi_{h,i}}=1|-)=1/\\{1+\exp(-2\xi_{h,i})\\}$; 2. Sample
$[\xi_{k,i}|-]\sim N(Q_{\xi_{k,i}}^{-1}\ell_{\xi_{k,i}},Q_{\xi_{k,i}}^{-1})$
where $Q_{\xi_{k,i}}=\eta_{k}^{2}/\sigma_{\epsilon}^{2}+1$ and
$\ell_{\xi_{k,i}}=\eta_{k}y_{k,i}^{F}/\sigma_{\epsilon}^{2}+m_{\xi_{k,i}}$; 3.
Sample $[\eta_{k}|-]\sim
N(Q_{\eta_{k}}^{-1}\ell_{\eta_{k}},Q_{\eta_{k}}^{-1})$ where
$Q_{\eta_{k}}=\sum_{i=1}^{n}\xi_{k,i}^{2}/\sigma_{\epsilon}^{2}+(\theta_{k}\sigma_{k}^{2})^{-1}$
and
$\ell_{\eta_{k}}=\sum_{i=1}^{n}\xi_{k,i}y_{k,i}^{F}/\sigma_{\epsilon}^{2}$; 4.
Rescale $\eta_{k}\rightarrow(\sum_{i=1}^{n}|\xi_{k,i}|/n)\eta_{k}$ and
$\bm{\xi}_{k}\rightarrow(n/\sum_{i=1}^{n}|\xi_{k,i}|)\bm{\xi}_{k}$ and update
$\beta_{k,i}=\xi_{k,i}\eta_{k}$; 5. Sample
$[\sigma_{k}^{-2}\mid-]\sim\mbox{Gamma}\left\\{a_{1}+1/2,a_{2}+\eta_{k}^{2}/(2\theta_{k})\right\\}$;
6. Sample
$[\nu_{k}\mid-]\sim\mbox{Beta}(1+\sum_{h=1}^{K}\mathbb{I}\\{z_{h}=k\\},\kappa+\sum_{h=1}^{K}\mathbb{I}\\{z_{h}>k\\})$
for $k=1,\ldots,K-1$ and update $\pi_{k}$ and $\omega_{k}$ from (6); 7. Sample
$[\kappa\mid-]\sim\mbox{Gamma}\\{a_{\kappa}+K-1,b_{\kappa}-\sum_{k=1}^{K-1}\log(1-\nu_{k})\\}$;
8. Sample $[z_{k}\mid-]$ from
$\mathbb{P}(z_{k}=h\mid-)\propto\begin{cases}\omega_{h}t_{2a_{1}}(\eta_{k};0,\sqrt{v_{0}a_{2}/a_{1}})&h\leq
k\\\ \omega_{h}t_{2a_{1}}(\eta_{k};0,\sqrt{a_{2}/a_{1}})&h>k\end{cases}$ where
$t_{d}(x;m,s)$ is the density of the $t$-distribution evaluated at $x$ with
mean $m$, standard deviation $s$, and degrees of freedom $d$; 9. Update
$\theta_{k}=1$ if $z_{k}>k$ and $\theta_{k}=v_{0}$ if $z_{k}\leq k$.

Algorithm 1 MCMC sampling steps for the ordered spike-and-slab prior

##  5 Simulation study

A simulation study is performed to assess model performance for point and
interval estimation of $Y_{i}$, $\bm{y}_{i}$, and $\\{\bm{\alpha}_{i}\\}$ and
uncertainty quantification for the number of nonparametric terms. We focus on
the linear template $\mathcal{H}_{0}=\mbox{span}\\{1,\tau\\}$ and present
results for the Nelson and Siegel, (1987) template with unknown $\gamma$ in
the supplementary material.

Synthetic functional data with $n=100$ curves and $m=25$ equally-spaced
observation points in $[0,1]$ are generated as follows. The parametric and
nonparametric factors are simulated as
$\alpha_{\ell,i}^{*}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,1)$ and
$\beta_{k,i}^{*}\stackrel{{\scriptstyle indep}}{{\sim}}N(0,1/(k+1)^{2})$,
respectively, for $k=1,\ldots,K_{true}$ and $K_{true}\in\\{0,1,3,8\\}$;
results for $K_{true}=1$ are nearly identical to the $K_{true}=3$ case and are
omitted, while $K_{true}=8$ is presented in the supplement. By design, the
variability in the parametric factors outweighs the variability in the
nonparametric factors. The parametric basis matrix $\bm{G}_{\gamma}^{*}$ is
constructed by evaluating $\\{g_{\ell}\\}$ at each observation point and QR-
decomposing the resulting matrix as in Section 2. For the nonparametric
functions $f_{k}^{*}$, we use orthogonal polynomials of degree $k+1$, which
are orthogonal to the linear template; the Nelson-Siegel version applies an
additional orthogonalization step. The error-free latent functions are
$\bm{Y}_{i}^{*}=\bm{G}_{\gamma}^{*}\bm{\alpha}_{i}^{*}+\bm{F}^{*}\bm{\beta}_{i}^{*}$
and the functional observations are generated as
$\bm{y}_{i}=\bm{Y}_{i}^{*}+\sigma^{*}\bm{\epsilon}_{i}^{*}$ where
$\sigma^{*}=\mbox{sd}(\bm{Y}_{i}^{*})/\mbox{RSNR}$ for sample standard
deviation $\mbox{sd}(\cdot)$, root signal-to-noise ratio $\mbox{RSNR}=3$, and
$\bm{\epsilon}_{i}^{*}\stackrel{{\scriptstyle
iid}}{{\sim}}N(\bm{0},\bm{I}_{m})$. This process was repeated to create 100
synthetic datasets.

We focus primarily on the PFFM, PFFM+gp, and SFFM. For all models, we assume
the conditionally Gaussian likelihood (4) with the hierarchical priors

$[\alpha_{\ell,i}\mid\sigma_{\alpha_{\ell}}]\stackrel{{\scriptstyle
indep}}{{\sim}}N(0,\sigma_{\alpha_{\ell}}^{2}),\quad\sigma_{\alpha_{\ell}}\stackrel{{\scriptstyle
iid}}{{\sim}}C^{+}(0,1),\quad p(\sigma_{\epsilon}^{2})\propto
1/\sigma_{\epsilon}^{2},$ (10)

and the SFFM hyperparameters $a_{1}=5$, $a_{2}=25$, and $v_{0}=0.001$, with an
upper bound $K=10$ on the number of nonparametric factors. Sensitivity
analyses were conducted for $(a_{1},a_{2})\in\\{(5,25),(5,50),(10,30)\\}$ and
$v_{0}\in\\{0.01,0.005,0.00025\\}$. The results for point and interval
predictions and estimates of $\bm{Y}_{i}^{*}$ and $\\{\alpha_{\ell,i}^{*}\\}$
are highly robust to these hyperparameters. Inference on $K^{*}$ is typically
robust with the exception of $v_{0}=0.00025$, for which the posterior of
$K^{*}$ becomes more sensitive to $(a_{1},a_{2})$. Hence, we select a larger
value of $v_{0}$.

First, we evaluate point prediction of the latent curves $\bm{Y}_{i}^{*}$ and
point estimation of the parametric factors $\\{\alpha_{\ell,i}^{*}\\}$.
Posterior expectations from each model are assessed using root mean squared
(prediction) error, and presented in Figure 4. Although the methods perform
similarly when the parametric model is true ($K_{true}=0$), the SFFM offers
substantial improvements in the semiparametric setting ($K_{true}>0$).
Notably, when $K_{true}>0$, the PFFM+gp offers better prediction than the
PFFM—as expected—but conversely cannot estimate the parametric factors
accurately.

Figure 4: Root mean squared (prediction) errors for $\bm{Y}^{*}$ (top) and
$\\{\alpha_{\ell,i}^{*}\\}$ (bottom) for $K_{true}=0$ (left) and $K_{true}=3$
(right). All methods perform similarly for $K_{true}=0$, while the SFFM
outperforms all competitors when $K_{true}>0$. In this setting, the PFFM+gp
offers better prediction than the PFFM, but cannot accurately estimate the
parametric factors.

Next, we assess prediction and credible intervals for $\bm{y}_{i}^{*}$ and
$\\{\alpha_{\ell,i}^{*}\\}$, respectively, where $\bm{y}_{i}^{*}$ is
distributed identically to $\bm{y}_{i}$. The 95% interval estimates are
evaluated using mean interval width and empirical coverage, which are
presented in Figure 5. Most notably, the SFFM provides uniformly narrower
interval estimates—and therefore more precise uncertainty quantification—with
approximately the correct nominal coverage. The PFFM prediction and credible
interval estimates are excessively wide when $K_{true}>0$, which is reflected
in Figure 1 and subsequently in Figure 7. The PFFM+gp interval estimates for
$\\{\alpha_{\ell,i}^{*}\\}$ are quite poor, especially for $K_{true}>0$, where
the interval estimates are excessively wide yet still fail to contain the true
parameters $\\{\alpha_{\ell,i}^{*}\\}$ at close to the nominal level. In
conjunction with Figure 2, these results show that functional confounding
adversely impacts both point estimation and uncertainty quantification.

Figure 5: Mean 95% prediction (top) and credible (bottom) interval widths for
$\bm{y}_{i}^{*}$ and $\\{\alpha_{\ell,i}^{*}\\}$, respectively, with empirical
coverage (annotations) for $K_{true}=0$ (left) and $K_{true}=3$ (right). The
SFFM provides narrow interval estimates with (nearly) the correct nominal
coverage in all cases. The PFFM intervals are excessively wide for
$K_{true}>0$, while the PFFM+gp performs poorly for
$\\{\alpha_{\ell,i}^{*}\\}$ in all cases.

Lastly, we study the posterior distribution of
$K^{*}=\sum_{k=1}^{K}\mathbb{I}\\{z_{k}>k\\}$ for estimating $K_{true}$ and
quantifying the necessity of the nonparametric terms. Since the PFFM and
PFFM+gp are not competitors for rank selection, we include a finite mixture
model variant of the SFFM (SFFM-fmm) that instead assumes
$\omega_{1:K}\sim\mbox{Dirichlet}(\kappa/K,\ldots,\kappa/K)$. We use the same
choice of $K$ as in the SFFM and set the prior expected number of factors to
be $\kappa=1$ to encourage fewer factors. Figure 6 compares the quality of the
posterior distribution $\mathbb{P}(K^{*}\mid\bm{y})$ for the SFFM and the
finite mixture alternative. Notably, the ordered spike-and-slab prior results
in much larger probabilities on the true rank $K_{true}$. Despite fixing
$\kappa=1$, the SFFM-fmm repeatedly overestimates the ranks: the simulation
average of $\mathbb{P}(K^{*}>K_{true}\mid\bm{y})$ is 0.30 for $K_{true}=0$,
0.23 for $K_{true}=3$, and 0.22 for $K_{true}=8$, while the comparable values
for the SFFM are 0.15, 0.16, and 0.15, respectively. Hence, the finite mixture
alternative introduces unnecessary parameters and relays deceptively strong
evidence for additional nonparametric terms.

Figure 6: Probability score $\mathbb{P}(K^{*}=K_{true}\mid\bm{y})$ for
$K_{true}=0$ (left) and $K_{true}=3$ (right). The proposed ordered spike-and-
slab prior provides better rank selection and inference than the finite
mixture alternative.

Additional simulation results in the supplement include $K_{true}=8$, the
Nelson and Siegel, (1987) template, and further assessments of nonparametric
detection via $\mathbb{P}(K^{*}>0\mid\bm{y})$, including ranked probability
scores for the distribution $\mathbb{P}(K^{*}\mid\bm{y})$. These results
confirm the results presented here.

##  6 Applications

###  6.1 Pinch force data

Human motor control is a critical area of research with implications for human
physiology, monitoring and mitigating muscle degeneration, and designing
robotic devices. Motor control data are recorded at high resolutions and often
modeled as functional data (Ramsay,, 2000; Goldsmith and Kitago,, 2016). Human
physiology and the laws of motion dictate a parametric template, which is
crucial for understanding the underlying process. However, these parametric
models may not adequately describe the observed data, which undermines the
interpretability of the key parameters.

We analyze human pinching data from Ramsay et al., (1995), which reports the
force measured by pinching the thumb and forefinger on opposing sides of a 6cm
force meter. The subject was instructed to maintain a background level of
constant force, then increase the pinching force to a predetermined maximum
level, and finally return to the original background level of constant force.
We use data from the fda package in R, which consists of $n=20$ replicate
force curves over time each with $m=151$ observations selected such that the
maximum of each curve occurred at 0.076 seconds. An example curve is in Figure
1.

To model the pinch force over time $\tau$, we adapt a parametric model from
Ramsay et al., (1995):

$g_{1}(\tau;\gamma_{i})=1,\quad g_{2}(\tau;\gamma_{i})=\exp[-(\log\tau-
c_{i})^{2}/\\{2\exp(\gamma_{i})\\}],$ (11)

where $\exp(c_{i})$ is the time of the maximum force and
$\gamma_{i}\in\mathbb{R}$ is a shape parameter. Ramsay et al., (1995) argue
that the unnormalized log-normal density $g_{2}$ matches the shape of the
observed data and offers plausible scientific explanations. For computational
convenience, we estimate $c_{i}$ as in Ramsay et al., (1995) by fitting a
quadratic regression in $\log(\tau)$ for the response variable
$\log(\bm{y}_{i})$ restricted to observations $y_{i}(\tau)>0.5$. The estimate
of each $c_{i}$ can be recovered from the estimated regression coefficients
and is subsequently treated as fixed.

We fit the PFFM, PFFM+gp, and SFFM to the data using the template (11). For
partial pooling among subjects, we specify a hierarchical prior on the shape
parameters: $\gamma_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}N(\mu_{\gamma},\sigma_{\gamma}^{2})$, $\mu_{\gamma}\sim N(0,10)$,
and $\sigma_{\gamma}\sim C^{+}(0,1)$, with the priors from (10) on the
remaining parameters. Since the nonlinear parameters $\gamma_{i}$ are curve-
specific, we construct each $\bm{G}_{\gamma_{i}}$ for $i=1,\ldots,n$ using a
QR decomposition as in Section 2 and modify the orthogonality constraint to be
$\bm{\bar{G}}^{\prime}\bm{f}_{k}=\bm{0}_{L}$ for
$\bm{\bar{G}}=n^{-1}\sum_{i=1}^{n}\bm{G}_{\gamma_{i}}$. This constraint no
longer preserves the posterior factorization of Corollaries 1-3, but
nonetheless maintains sufficient distinctness between the parametric and
nonparametric terms. The SFFM hyperparameters are set to $a_{1}=5$,
$a_{2}=25$, $v_{0}=0.001$, and $K=10$ as in Section 5.

An example of the fitted values with 95% simultaneous prediction bands for the
PFFM and SFFM is in Figure 1. Although the PFFM captures the general shape of
the data, it suffers from clear bias around the peak and produces
unnecessarily wide prediction bands. The SFFM corrects both issues: the bias
is removed and the prediction bands are more precise. Notably, the fitted SFFM
curve preserves the same general shape as the PFFM and avoids overfitting
despite the increase in modeling complexity.

Posterior uncertainty quantification is also more precise under the SFFM for
the parametric factors $\\{\alpha_{\ell,i}\\}$. Figure 7 presents the
posterior standard deviations of $\\{\alpha_{\ell,i}\\}_{i=1}^{n}$ for each
$\ell=1,\ldots,L$ under the PFFM and SFFM. Although the posterior expectations
are similar (see Figure 2), the posterior standard deviations are
substantially smaller under the SFFM. This reduction is even more dramatic
compared to the PFFM+gp, which is omitted from Figure 7 due to excessively
large posterior standard deviations that range from 1.42 to 4.43. For model
(11), the linear coefficients $\\{\alpha_{\ell,i}\\}$ determine the maximum of
the force curve, which is the most prominent feature in the data. The SFFM
provides more precise posterior inference for these key parametric factors.

Figure 7: Posterior standard deviations of $\\{\alpha_{\ell,i}\\}_{i=1}^{n}$
for $\ell=1,\ldots,L$ in the parametric (PFFM) and semiparametric functional
factor model (SFFM) for the pinch force data (left) and yield curve data
(right). There is a clear and consistent reduction in posterior uncertainty
from the PFFM to the SFFM for these key parametric factors.

Additional evidence in favor of the SFFM is presented in Table 2, which
estimates the posterior distribution of the effective number of nonparametric
terms $K^{*}$. Clearly, there is substantially posterior probability for at
least two nonparametric terms, with limited evidence that more than three
terms are needed. WAIC similarly favors the SFFM: the estimates are $-922$,
$-2693$, and $-3005$ for the PFFM, PFFM+gp, and SFFM, respectively.

$k$ | 0 | 1 | 2 | 3 | 4 | 5 | $\geq 6$
---|---|---|---|---|---|---|---
$\mathbb{P}(K^{*}=k\mid\bm{y})$ | 0.003 | 0.147 | 0.707 | 0.135 | 0.007 | 0.001 | 0
Table 2: Posterior probabilities $\mathbb{P}(K^{*}=k\mid\bm{y})$ for the pinch
force data. There is strong evidence for a nonparametric component.

We highlight the excellent MCMC performance of the SFFM. For all models, we
retain 10000 MCMC samples after discarding 5000 iterations as a burn-in.
Traceplots of $\\{Y_{i}\\}$, $\\{\alpha_{\ell,i}\\}$, and $K^{*}$ (not shown)
demonstrate excellent mixing and suggest convergence. The computing time per
1000 iterations is 1.59s, 4.35s, and 1.98s for the PFFM, PFFM+gp, and SFFM,
respectively (using R on a MacBook Pro, 2.8 GHz Intel Core i7). Comparing the
effective sample sizes for all $\\{\alpha_{\ell,i}\\}$, the first quartiles
are 10000, 1361, and 9508 for the PFFM, PFFM+gp, and SFFM, respectively. In
aggregate, the SFFM not only provides superior modeling capabilities, but also
rivals the simpler PFFM in computing time and MCMC efficiency—while the
PFFM+gp lags far behind in both.

###  6.2 Dynamic yield curves

Yield curves play a fundamental role in economic and financial analyses: they
provide essential information about current and future economic conditions,
including inflation, business cycles, and monetary policies, and are used to
price fixed-income securities and construct forward curves. The yield curve
$Y_{i}(\tau)$ describes how interest rates vary as a function of the length of
the borrowing period, or time to maturity $\tau$, at each time $i$. Naturally,
yield curves can be modeled as functional data that evolve dynamically over
time.

Yield curve models most commonly employ the parametric Nelson and Siegel,
(1987) basis:

$g_{1}(\tau;\gamma)=1,\quad
g_{2}(\tau;\gamma)=\\{1-\exp(-\tau\gamma)\\}/(\tau\gamma),\quad
g_{3}(\tau;\gamma)=g_{2}(\tau;\gamma)-\exp(-\tau\gamma)$ (12)

where $g_{1}$ is the _level_ , $g_{2}$ is the _slope_ , and $g_{3}$ is the
_curvature_. The nonlinear parameter $\gamma>0$ is commonly treated as fixed,
such as $\gamma=0.0609$ (Diebold and Li,, 2006), but otherwise may be
estimated. We use a weakly informative Gamma prior for $\gamma$ with prior
mean 0.0609 and prior variance 0.5. Both the PFFM and the SFFM with this prior
on $\gamma$ are favored by WAIC over their respective counterparts with fixed
$\gamma=0.0609$.

To capture the yield curve dynamics, we model the parametric Nelson-Siegel
factors as an AR(1):

$\alpha_{\ell,i}=\mu_{\ell}+\phi_{\ell}(\alpha_{\ell,i-1}-\mu_{\ell})+\zeta_{\ell,i},\quad\zeta_{\ell,i}\stackrel{{\scriptstyle
indep}}{{\sim}}N(0,\sigma_{\zeta_{\ell}}^{2}).$ (13)

The PFFM with (12) and (13) is also known as the _dynamic Nelson-Siegel model_
(Diebold and Li,, 2006). The dynamic Nelson-Siegel factors
$\\{\alpha_{\ell,i}\\}$ may be viewed as the state variables in a dynamic
linear model. Using the orthogonality constraints on $\\{g_{\ell}\\}$ and
$\\{f_{k}\\}$, Corollary 2 provides a convenient and low-dimensional state
simulation algorithm for efficient joint sampling of the dynamic factors
$\\{\alpha_{\ell,i}\\}$; see the supplementary material for details. Similar
simplifications are _not_ available for the PFFM+gp, which is omitted.

The dynamic model (13) is accompanied by a diffuse prior
$\mu_{\ell}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,10^{6})$ and a weakly
informative prior $(\phi_{\ell}+1)/2\stackrel{{\scriptstyle
iid}}{{\sim}}\mbox{Beta}(5,2)$—which ensures stationarity of the dynamic
factors $\\{\alpha_{\ell,i}\\}$ and therefore $\\{Y_{i}\\}$ and
$\\{\bm{y}_{i}\\}$—along with $\sigma_{\zeta_{\ell}}\stackrel{{\scriptstyle
iid}}{{\sim}}C^{+}(0,1)$ for $\ell=1,\ldots,L$. In addition, we generalize (2)
to accommodate stochastic volatility in the observation error variance,
$\epsilon_{i,j}\stackrel{{\scriptstyle
indep}}{{\sim}}N(0,\sigma_{\epsilon_{i}}^{2})$, with
$\log\sigma_{\epsilon_{i}}^{2}\sim\mbox{AR(1)}$, which is an essential
component in many economic and financial models (Kim et al.,, 1998). For these
AR(1) parameters, we adopt the priors and sampling algorithm from Kowal,
(2020).

We evaluate the suitability of the Nelson-Siegel basis for modeling monthly
unsmoothed Fama-Bliss US government bond yields from 2000-2009 ($n=120$)
provided by Van Dijk et al., (2014). These data are available for maturities
of 3, 6, 9, 12, 15, 18, 21, 24, 30, 36, 48, 60, 72, 84, 96, 108 and 120 months
($m=17$). The inferential targets are the latent curves $Y_{i}$, the dynamic
Nelson-Siegel factors $\\{\alpha_{i,\ell}\\}$, and the effective number of
nonparametric terms $K^{*}$. The SFFM hyperparameter and MCMC specifications
from Section 6.1 are adopted here, again with excellent mixing and convergence
for these key parameters.

Figure 1 shows the PFFM and SFFM fitted values with 95% simultaneous
prediction bands for the yield curve in September 2008 during the onset of the
Great Recession. The PFFM produces a reasonable shape for the yield curve,
yet—as with the pinch force data—suffers from clear bias and overconservative
prediction bands. The SFFM corrects these deficiencies without distorting the
general shape of the curve or overfitting to the data.

Similar results are obtained in Figure 7 for the dynamic Nelson-Siegel
factors. The SFFM offers substantial reductions in posterior standard
deviation for all three dynamic Nelson-Siegel factors. These factors are
fundamental for describing the shape of the yield curve; reducing the
posterior uncertainty is a crucial advantage of the SFFM. Importantly, the
simulation study confirms that the reduced posterior uncertainty
quantification from the SFFM nonetheless retains valid calibration, or more
specifically, correct nominal coverage of the posterior credible intervals.

We summarize the posterior distribution $\mathbb{P}(K^{*}\mid\bm{y})$ of the
effective number of nonparametric terms in Table 3. The evidence for the
nonparametric factors is moderate: we estimate
$\mathbb{P}(K^{*}>0\mid\bm{y})=0.09$.

$k$ | 0 | 1 | 2 | $\geq 3$
---|---|---|---|---
$\mathbb{P}(K^{*}=k|\bm{y})$ | 0.910 | 0.087 | 0.003 | 0
Table 3: Posterior probabilities $\mathbb{P}(K^{*}=k\mid\bm{y})$ for the yield
curve data. There is moderate evidence for a nonparametric component.

There are other compelling reasons to include the nonparametric factors in a
dynamic Nelson-Siegel model. First, the uncertainty quantification for
$\\{Y_{i}\\}$ and $\\{\alpha_{\ell,i}\\}$ is more precise in the SFFM than in
the PFFM (see Figures 1 and 7). Second, WAIC decisively prefers the SFFM
(-10229) over the PFFM (-5516), which suggests potential improvements in out-
of-sample predictive capabilities. Third, the additional computational burdens
of the SFFM are minimal. Due to the orthogonality constraints and Corollory 2,
the sampling steps for $\\{\alpha_{\ell,i}\\}$ and associated parameters are
identical for the PFFM and the SFFM, which permits full and robust model
development for the parametric factors. Further, the computing times per 1000
iterations are similar: 21s and 23s for the PFFM and SFFM, respectively (using
R on a MacBook Pro, 2.8 GHz Intel Core i7).

##  7 Discussion

We proposed a Bayesian semiparametric model for functional data. The
semiparametric functional factor model (SFFM) augmented a parametric template
with an infinite-dimensional nonparametric functional basis. The nonparametric
basis was treated as unknown and learned from the data to correct the biases
of the parametric template while appropriately incorporating relevant
uncertainties into the posterior distribution. Distinctness between the
parametric and nonparametric terms was achieved by conditioning upon an
orthogonality constraint, which simultaneously prevented _functional
confounding_ and admitted highly convenient simplifications for efficient MCMC
sampling. The nonparametric component was regularized with an ordered spike-
and-slab prior, which implicitly provided rank selection for infinite-
dimensional models and satisfied appealing theoretical properties. This prior
was accompanied by a parameter expansion scheme customized to boost MCMC
efficiency, and is broadly applicable for Bayesian factor models. Our analyses
of synthetic data, human motor control data, and dynamic interest rates
demonstrated clear advantages of the semiparametric modeling framework
relative to both parametric and Gaussian process alternatives. The proposed
approach eliminated bias, reduced excessive posterior and predictive
uncertainty, and provided reliable inference on the effective number of
nonparametric terms—all with minimal computational costs.

There are several promising extensions that remain for future work. First,
formulation of the parametric template in (1) can be generalized, for example
to include a functional regression term with scalar or functional covariates.
Second, the ordered spike-and-slab prior currently uses independent and
identically distributed variables in the parameter expansion (9). Adaptations
to include dependence among these variables, such as regression models,
clustering, and spatio-temporal dependence, would broaden the applicability of
the prior. Lastly, we studied only a small subset of many possible parametric
templates. For applications that rely on such parametric models, the proposed
semiparametric modeling framework can directly assess the adequacy of these
models—and perhaps suggest improvements.

### Acknowledgements

Research was sponsored by the Army Research Office and was accomplished under
Grant Number W911NF-20-1-0184 (Kowal). The views and conclusions contained in
this document are those of the authors and should not be interpreted as
representing the official policies, either expressed or implied, of the Army
Research Office or the U.S. Government. The U.S. Government is authorized to
reproduce and distribute reprints for Government purposes notwithstanding any
copyright notation herein.

## References

  * Bhattacharya and Dunson, (2011) Bhattacharya, A. and Dunson, D. B. (2011).  Sparse Bayesian infinite factor models.  Biometrika, 98(2):291–306.
  * Bianchi et al., (2009) Bianchi, F., Mumtaz, H., and Surico, P. (2009).  The great moderation of the term structure of UK interest rates.  Journal of Monetary Economics, 56(6):856–871.
  * Canale et al., (2017) Canale, A., Lijoi, A., Nipoti, B., and Prünster, I. (2017).  On the Pitman-Yor process with spike and slab base measure.  Biometrika, 104(3):681–697.
  * Castillo and van der Vaart, (2012) Castillo, I. and van der Vaart, A. (2012).  Needles and straw in a haystack: Posterior concentration for possibly sparse sequences.  The Annals of Statistics, 40(4):2069–2101.
  * Chen et al., (2019) Chen, S. T., Xiao, L., and Staicu, A.-M. (2019).  Model Testing for Generalized Scalar-on-Function Linear Models.  arXiv preprint arXiv:1906.04889.
  * Diebold and Li, (2006) Diebold, F. X. and Li, C. (2006).  Forecasting the term structure of government bond yields.  Journal of Econometrics, 130(2):337–364.
  * Durbin and Koopman, (2002) Durbin, J. and Koopman, S. J. (2002).  A simple and efficient simulation smoother for state space time series analysis.  Biometrika, 89(3):603–615.
  * Goldsmith and Kitago, (2016) Goldsmith, J. and Kitago, T. (2016).  Assessing systematic effects of stroke on motor control by using hierarchical function-on-scalar regression.  Journal of the Royal Statistical Society. Series C: Applied Statistics, 65(2):215–236.
  * Griffiths and Ghahramani, (2011) Griffiths, T. L. and Ghahramani, Z. (2011).  The Indian Buffet Process: An Introduction and Review.  Journal of Machine Learning Research, 12(4).
  * Heckman and Ramsay, (2000) Heckman, N. E. and Ramsay, J. O. (2000).  Penalized regression with model-based penalties.  Canadian Journal of Statistics, 28(2):241–258.
  * Ishwaran and Rao, (2005) Ishwaran, H. and Rao, J. S. (2005).  Spike and slab variable selection: Frequentist and Bayesian strategies.  Annals of Statistics, 33(2):730–773.
  * Kim et al., (1998) Kim, S., Shephard, N., and Chib, S. (1998).  Stochastic Volatility: Likelihood Inference and Comparison with ARCH Models.  Review of Economic Studies, 65(3):361–393.
  * Kowal, (2020) Kowal, D. R. (2020).  Dynamic Regression Models for Time-Ordered Functional Data.  Bayesian Analysis.
  * Kowal and Bourgeois, (2020) Kowal, D. R. and Bourgeois, D. C. (2020).  Bayesian Function-on-Scalars Regression for High-Dimensional Data.  Journal of Computational and Graphical Statistics, 29(3):1–10.
  * Lee et al., (2018) Lee, W., Miranda, M. F., Rausch, P., Baladandayuthapani, V., Fazio, M., Downs, J. C., and Morris, J. S. (2018).  Bayesian semiparametric functional mixed models for serially correlated functional data, with application to glaucoma data.  Journal of the American Statistical Association.
  * Legramanti et al., (2020) Legramanti, S., Durante, D., and Dunson, D. B. (2020).  Bayesian cumulative shrinkage for infinite factorizations.  Biometrika, 107(3):745–752.
  * Mikulich et al., (2003) Mikulich, S. K., Zerbe, G. O., Jones, R. H., and Crowley, T. J. (2003).  Comparing linear and nonlinear mixed model approaches to cosinor analysis.  Statistics in Medicine, 22(20):3195–3211.
  * Molenberghs and Verbeke, (2000) Molenberghs, G. and Verbeke, G. (2000).  Linear mixed models for longitudinal data.  Springer.
  * Neal, (2003) Neal, R. M. (2003).  Slice sampling.  Annals of Statistics, pages 705–741.
  * Nelson and Siegel, (1987) Nelson, C. R. and Siegel, A. F. (1987).  Parsimonious Modeling of Yield Curves.  The Journal of Business, 60(4):473.
  * Ohn and Kim, (2021) Ohn, I. and Kim, Y. (2021).  Posterior Consistency of Factor Dimensionality in High-Dimensional Sparse Factor Models.  Bayesian Analysis, 1(1):1–24.
  * Ramsay, (2000) Ramsay, J. O. (2000).  Functional components of variation in handwriting.  Journal of the American Statistical Association, 95(449):9–15.
  * Ramsay and Dalzell, (1991) Ramsay, J. O. and Dalzell, C. J. (1991).  Some tools for functional data analysis.  Journal of the Royal Statistical Society: Series B (Methodological), 53(3):539–561.
  * Ramsay et al., (1995) Ramsay, J. O., Wang, X., and Flanagan, R. (1995).  A functional data analysis of the pinch force of human fingers.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 44(1):17–30.
  * Reich et al., (2006) Reich, B. J., Hodges, J. S., and Zadnik, V. (2006).  Effects of residual smoothing on the posterior of the fixed effects in disease‐mapping models.  Biometrics, 62(4):1197–1206.
  * Rockova, (2018) Rockova, V. (2018).  Bayesian estimation of sparse signals with a continuous spike and slab prior.  The Annals of Statistics, 46(1):401–437.
  * Sang et al., (2017) Sang, P., Wang, L., and Cao, J. (2017).  Parametric functional principal component analysis.  Biometrics, 73(3):802–810.
  * Scarpa and Dunson, (2009) Scarpa, B. and Dunson, D. B. (2009).  Bayesian hierarchical functional data analysis via contaminated informative priors.  Biometrics, 65(3):772–780.
  * Scarpa and Dunson, (2014) Scarpa, B. and Dunson, D. B. (2014).  Enriched stick-breaking processes for functional data.  Journal of the American Statistical Association, 109(506):647–660.
  * Scheipl et al., (2012) Scheipl, F., Fahrmeir, L., and Kneib, T. (2012).  Spike-and-slab priors for function selection in structured additive regression models.  Journal of the American Statistical Association, 107(500):1518–1532.
  * Teh et al., (2007) Teh, Y. W., Grür, D., and Ghahramani, Z. (2007).  Stick-breaking construction for the Indian buffet process.  In Artificial Intelligence and Statistics, pages 556–563. PMLR.
  * Van Dijk et al., (2014) Van Dijk, D., Koopman, S. J., der Wel, M., and Wright, J. H. (2014).  Forecasting interest rates with shifting endpoints.  Journal of Applied Econometrics, 29(5):693–712.
  * Wagner et al., (2020) Wagner, A. B., Hill, E. L., Ryan, S. E., Sun, Z., Deng, G., Bhadane, S., Martinez, V. H., Wu, P., Li, D., Anand, A., Acharya, J., and Matteson, D. S. (2020).  Social Distancing Has Merely Stabilized COVID-19 in the US.  Stat, page e302.
  * Wand et al., (2011) Wand, M. P., Ormerod, J. T., Padoan, S. A., and Frührwirth, R. (2011).  Mean field variational Bayes for elaborate distributions.  Bayesian Analysis, 6(4):847–900.
  * Welham et al., (2006) Welham, S. J., Cullis, B. R., Kenward, M. G., and Thompson, R. (2006).  The analysis of longitudinal data using mixed model L‐splines.  Biometrics, 62(2):392–401.

##  8 Proofs

######  Lemma 4.

When $\bm{G}_{\gamma}^{\prime}\bm{f}_{k}=\bm{0}_{L}$ for all $k=1,2,\ldots$
and conditional on $\sigma_{\epsilon}^{2}$, the likelihood (4) factorizes:
$p(\bm{y}\mid\gamma,\\{\bm{\alpha}_{i}\\},\bm{F},\\{\bm{\beta}_{i}\\})=p_{0}(\bm{y}\mid\gamma,\\{\bm{\alpha}_{i}\\})p_{1}(\bm{y}\mid\bm{F},\\{\bm{\beta}_{i}\\})$,
where $p_{0}$ depends _only_ on the parametric terms $\gamma$ and
$\\{\bm{\alpha}_{i}\\}$ and $p_{1}$ depends _only_ on the nonparametric terms
$\bm{F}$ and $\\{\bm{\beta}_{i}\\}$.

###### Proof.

Observing that
$\|\bm{y}_{i}-\bm{G}_{\gamma}\bm{\alpha}_{i}-\bm{F}\bm{\beta}_{i}\|^{2}=\|\bm{y}_{i}\|^{2}+\|\bm{G}_{\gamma}\bm{\alpha}_{i}\|^{2}+\|\bm{F}\bm{\beta}_{i}\|^{2}-2\bm{\alpha}_{i}\bm{G}_{\gamma}^{\prime}\bm{y}_{i}-2\bm{\alpha}_{i}\bm{G}_{\gamma}^{\prime}\bm{F}\bm{\beta}_{i}-2\bm{\beta}_{i}\bm{F}^{\prime}\bm{y}_{i}=\|\bm{y}_{i}\|^{2}+\big{(}\|\bm{G}_{\gamma}\bm{\alpha}_{i}\|^{2}-2\bm{\alpha}_{i}\bm{G}_{\gamma}^{\prime}\bm{y}_{i}\big{)}+\big{(}\|\bm{F}\bm{\beta}_{i}\|^{2}-2\bm{\beta}_{i}\bm{F}^{\prime}\bm{y}_{i}\big{)}$,
the joint likelihood factorizes into

$\displaystyle p(\bm{y}|-)$
$\displaystyle=\prod_{i=1}^{n}p(\bm{y}_{i}|-)=\prod_{i=1}^{n}(2\pi\sigma_{\epsilon}^{2})^{-m/2}\exp\\{-\|\bm{y}_{i}-\bm{G}_{\gamma}\bm{\alpha}_{i}-\bm{F}\bm{\beta}_{i}\|^{2}/(2\sigma_{\epsilon}^{2})\\}$
$\displaystyle=\prod_{i=1}^{n}(2\pi\sigma_{\epsilon}^{2})^{-m/2}\exp[-\\{\|\bm{y}_{i}\|^{2}+\big{(}\|\bm{G}_{\gamma}\bm{\alpha}_{i}\|^{2}-2\bm{\alpha}_{i}\bm{G}_{\gamma}^{\prime}\bm{y}_{i}\big{)}+\big{(}\|\bm{F}\bm{\beta}_{i}\|^{2}-2\bm{\beta}_{i}\bm{F}^{\prime}\bm{y}_{i}\big{)}\\}/(2\sigma_{\epsilon}^{2})]$
$\displaystyle=(2\pi\sigma_{\epsilon}^{2})^{-nm/2}\exp\Big{\\{}-\sum_{i=1}^{n}\|\bm{y}_{i}\|^{2}/(2\sigma_{\epsilon}^{2})\Big{\\}}\exp\Big{\\{}-\sum_{i=1}^{n}\big{(}\|\bm{G}_{\gamma}\bm{\alpha}_{i}\|^{2}-2\bm{\alpha}_{i}^{\prime}\bm{G}_{\gamma}^{\prime}\bm{y}_{i}\big{)}/(2\sigma_{\epsilon}^{2})\Big{\\}}$
$\displaystyle\quad\quad\quad\exp\Big{\\{}-\sum_{i=1}^{n}\big{(}\|\bm{F}\bm{\beta}_{i}\|^{2}-2\bm{\beta}_{i}^{\prime}\bm{F}^{\prime}\bm{y}_{i}\big{)}/(2\sigma_{\epsilon}^{2})\Big{\\}}$

which produces the necessary product factorization. ∎

######  Corollary 5.

Under model (4) with prior independence
$p(\\{\bm{\alpha}_{i}\\},\\{\bm{\beta}_{i}\\})=p(\\{\bm{\alpha}_{i}\\})p(\\{\bm{\beta}_{i}\\})$
and $\bm{G}_{\gamma}^{\prime}\bm{f}_{k}=\bm{0}_{L}$ for all $k=1,2,\ldots$,
the parametric and nonparametric factors are _a posteriori_ independent:
$p(\\{\bm{\alpha}_{i}\\},\\{\bm{\beta}_{i}\\}|\bm{y},\sigma_{\epsilon}^{2},\bm{F},\gamma)=p(\\{\bm{\alpha}_{i}\\}|\bm{y},\sigma_{\epsilon}^{2},\bm{F},\gamma)p(\\{\bm{\beta}_{i}\\}|\bm{y},\sigma_{\epsilon}^{2},\bm{F},\gamma)$.

###### Proof.

The result follows as a direct consequence of (i) the prior independence
assumption and (ii) the likelihood factorizations from Lemma 1. ∎

######  Lemma 5.

Suppose $\bm{G}_{\gamma}$ is $L\times m$ with rank $L$. Denote
$\bm{\psi}_{k}^{0}\sim
N\left(\bm{Q}_{\psi_{k}}^{-1}\bm{\ell}_{\psi_{k}},\bm{Q}_{\psi_{k}}^{-1}\right)$
with $\bm{f}_{k}^{0}=\bm{B}\bm{\psi}_{k}^{0}$. The shifted term
$\bm{f}_{k}=\bm{B}\bm{\psi}_{k}$ with
$\bm{\psi}_{k}=\bm{\psi}_{k}^{0}-\bm{Q}_{\psi_{k}}^{-1}\bm{B}^{\prime}\bm{G}_{\gamma}\left(\bm{G}_{\gamma}^{\prime}\bm{B}\bm{Q}_{\psi_{k}}^{-1}\bm{B}^{\prime}\bm{G}_{\gamma}\right)^{-1}\bm{G}_{\gamma}^{\prime}\bm{B}\bm{\psi}_{k}^{0}$
satisfies $\bm{f}_{k}\stackrel{{\scriptstyle
d}}{{=}}[\bm{f}_{k}^{0}\mid\bm{G}_{\gamma}^{\prime}\bm{f}_{k}^{0}=\bm{0}]$ and
$\mathbb{P}(\bm{G}_{\gamma}^{\prime}\bm{f}_{k}=\bm{0})=1$.

###### Proof.

To show the result, we apply a modification of Theorem 1 in Kowal, (2020). The
full conditional distribution of the unconstrained term is $\bm{f}_{k}^{0}\sim
N\big{(}\bm{\mu}_{f_{k}^{0}},\bm{\Sigma}_{f_{k}^{0}}\big{)}$ with
$\bm{\mu}_{f_{k}^{0}}=\bm{B}\bm{Q}_{\psi_{k}}^{-1}\bm{\ell}_{\psi_{k}}$ and
$\bm{\Sigma}_{f_{k}^{0}}=\bm{B}\bm{Q}_{\psi_{k}}^{-1}\bm{B}^{\prime}$. Since
$\bm{f}_{k}^{0}$ is Gaussian, so is $\bm{G}_{\gamma}^{\prime}\bm{f}_{k}^{0}$,
and consequently
$[\bm{f}_{k}^{0}|\bm{G}_{\gamma}^{\prime}\bm{f}_{k}^{0}=\bm{0}]\sim
N\big{(}\bm{\mu}_{f_{k}},\bm{\Sigma}_{f_{k}}\big{)}$ with
$\bm{\mu}_{f_{k}}=\bm{\mu}_{f_{k}^{0}}-\bm{\Sigma}_{f_{k}^{0}}\bm{G}_{\gamma}(\bm{G}_{\gamma}^{\prime}\bm{\Sigma}_{f_{k}^{0}}\bm{G}_{\gamma})^{-1}\bm{G}_{\gamma}^{\prime}\bm{\mu}_{f_{k}^{0}}$
and
$\bm{\Sigma}_{f_{k}}=\bm{\Sigma}_{f_{k}^{0}}-\bm{\Sigma}_{f_{k}^{0}}\bm{G}_{\gamma}(\bm{G}_{\gamma}^{\prime}\bm{\Sigma}_{f_{k}^{0}}\bm{G}_{\gamma})^{-1}\bm{\Sigma}_{f_{k}^{0}}$.
Direct calculation shows that the shifted term $\bm{f}_{k}$ has the same mean
and covariance. Pre-multiplying a random variable with distribution
$[\bm{f}_{k}^{0}|\bm{G}_{\gamma}^{\prime}\bm{f}_{k}^{0}=\bm{0}]$ by
$\bm{G}_{\gamma}^{\prime}$ produces a random variable with degenerate
distribution at $\bm{0}_{L}$. ∎

######  Corollary 6.

Under model (4) and subject to $\bm{G}_{\gamma}^{\prime}\bm{f}_{k}=\bm{0}_{L}$
for all $k$ and $\bm{G}_{\gamma}^{\prime}\bm{G}_{\gamma}=\bm{I}_{L}$, the
likelihood for the parametric factors is proportional to
$p(\bm{y}\mid\gamma,\\{\bm{\alpha}_{i}\\},\bm{F},\\{\bm{\beta}_{i}\\})\propto
p_{0,G}(\bm{y}\mid\gamma,\\{\bm{\alpha}_{i}\\})$ where $p_{0,G}$ is the
likelihood defined by
$[\bm{g}_{\ell;\gamma}^{\prime}\bm{y}_{i}\mid-]\stackrel{{\scriptstyle
indep}}{{\sim}}N(\alpha_{\ell,i},\sigma_{\epsilon}^{2})$ for
$\ell=1,\ldots,L$.

###### Proof.

From the proof of Lemma 1 and considering terms involving only the parametric
factors, we have

$\displaystyle p(\bm{y}|-)$
$\displaystyle\propto\exp\Big{\\{}-\sum_{i=1}^{n}\big{(}\|\bm{G}_{\gamma}\bm{\alpha}_{i}\|^{2}-2\bm{\alpha}_{i}\bm{G}_{\gamma}^{\prime}\bm{y}_{i}\big{)}/(2\sigma_{\epsilon}^{2})\Big{\\}}$
$\displaystyle=\exp\Big{\\{}-\sum_{i=1}^{n}\big{(}\|\bm{\alpha}_{i}\|^{2}-2\bm{\alpha}_{i}^{\prime}\bm{G}_{\gamma}^{\prime}\bm{y}_{i}\big{)}/(2\sigma_{\epsilon}^{2})\Big{\\}}$
$\displaystyle\propto p_{0,G}(\bm{y}\mid\gamma,\\{\bm{\alpha}_{i}\\})$

using the orthonormality $\bm{G}_{\gamma}^{\prime}\bm{G}_{\gamma}=\bm{I}_{L}$.
∎

######  Corollary 7.

Under model (4) and subject to $\bm{G}_{\gamma}^{\prime}\bm{f}_{k}=\bm{0}_{L}$
for all $k$ and $\bm{F}^{\prime}\bm{F}=\bm{I}_{\infty}$, the likelihood for
the nonparametric factors is proportional to
$p(\bm{y}\mid\gamma,\\{\bm{\alpha}_{i}\\},\bm{F},\\{\bm{\beta}_{i}\\})\propto
p_{1,F}(\bm{y}\mid\bm{F},\\{\bm{\beta}_{i}\\})$ where $p_{1,F}$ is the
likelihood defined by
$[\bm{f}_{k}^{\prime}\bm{y}_{i}\mid-]\stackrel{{\scriptstyle
indep}}{{\sim}}N(\beta_{k,i},\sigma_{\epsilon}^{2})$ for all $k$.

###### Proof.

The result follows using the same arguments as in the proof of Corollary 2. ∎

######  Proposition 5.

For $\varepsilon>0$ and fixed $\theta_{0}$, let
$\mathbb{B}_{\varepsilon}(\theta_{0})=\\{\theta_{k}:|\theta_{k}-\theta_{0}|<\varepsilon\\}$.
Prior (5)-(6) implies that
$\mathbb{P}(|\theta_{k}-\theta_{0}|\leq\varepsilon)<\mathbb{P}(|\theta_{k+1}-\theta_{0}|\leq\varepsilon)$
whenever
$P_{slab}\\{\mathbb{B}_{\varepsilon}(\theta_{0})\\}<P_{spike}\\{\mathbb{B}_{\varepsilon}(\theta_{0})\\}$.

###### Proof.

The proof proceeds directly:

$\displaystyle\mathbb{P}(|\theta_{k}-\theta_{0}|\leq\varepsilon)$
$\displaystyle=\mathbb{E}[P_{k}\\{\mathbb{B}_{\varepsilon}(\theta_{0})\\}]=\mathbb{E}[(1-\pi_{k})P_{slab}\\{\mathbb{B}_{\varepsilon}(\theta_{0})\\}+\pi_{k}P_{spike}\\{\mathbb{B}_{\varepsilon}(\theta_{0})\\}]$
$\displaystyle=\mathbb{E}\\{(1-\pi_{k})\\}P_{slab}\\{\mathbb{B}_{\varepsilon}(\theta_{0})\\}+\mathbb{E}(\pi_{k})P_{spike}\\{\mathbb{B}_{\varepsilon}(\theta_{0})\\}$
$\displaystyle=P_{slab}\\{\mathbb{B}_{\varepsilon}(\theta_{0})\\}\\{\kappa/(1+\kappa)\\}^{k}+P_{spike}\\{\mathbb{B}_{\varepsilon}(\theta_{0})\\}[1-\\{\kappa/(1+\kappa)\\}^{k}]$
$\displaystyle=P_{spike}\\{\mathbb{B}_{\varepsilon}(\theta_{0})\\}+\\{\kappa/(1+\kappa)\\}^{k}[P_{slab}\\{\mathbb{B}_{\varepsilon}(\theta_{0})\\}-P_{spike}\\{\mathbb{B}_{\varepsilon}(\theta_{0})\\}]$

noting
$\mathbb{E}(\pi_{k})=1-\\{\iota\kappa/(\iota+\iota\kappa)\\}^{k}=1-\\{\kappa/(1+\kappa)\\}^{k}$.
Since $\\{\kappa/(1+\kappa)\\}^{k}\in[0,1]$ is decreasing in $k$ and the
remaining terms are invariant to $k$, the result follows. ∎

######  Corollary 8.

For $\varepsilon>0$,
$\mathbb{P}(|\eta_{k}|\leq\varepsilon)<\mathbb{P}(|\eta_{k+1}|\leq\varepsilon)$
whenever $v_{0}<1$.

###### Proof.

Corollary 4 follows from Proposition 1 by observing that
$[\eta_{k}|\pi_{k}]\sim(1-\pi_{k})t_{2a_{1}}(0,\sqrt{a_{2}/a_{1}})+\pi_{k}t_{2a_{1}}(0,\sqrt{v_{0}a_{2}/a_{1}})$,
where $t_{d}(m,s)$ denotes a $t$-distribution with mean $m$, standard
deviation $s$, and degrees of freedom $d$. This mixture is a special case of
(5), and the densities of the $t$-distributions place greater mass near zero
when the scale parameter is smaller. The condition $v_{0}<1$ ensures that the
spike distribution is indeed more concentrated around zero. ∎

######  Proposition 6.

The CUSP (6) satisfies $(1-\pi_{k})=\mu_{(k)}$ where
$\mu_{(1)}>\cdots>\mu_{(K)}$ are the ordered (slab) probabilities from the Teh
et al., (2007) stick-breaking construction of the IBP, i.e.,
$\mu_{(k)}=\prod_{\ell=1}^{k}\nu_{\ell}^{\prime}$ with
$\nu_{\ell}^{\prime}\stackrel{{\scriptstyle
iid}}{{\sim}}\mbox{Beta}(\iota\kappa,\iota).$

###### Proof.

Recall that IBPs are obtained by establishing the conditional distributions
$[b_{j,k}\mid\mu_{k}]\stackrel{{\scriptstyle
indep}}{{\sim}}\mbox{Bernoulli}(\mu_{k})$ and
$[\mu_{k}]\stackrel{{\scriptstyle
iid}}{{\sim}}\mbox{Beta}(\iota\kappa/K,\iota)$ and then integrating over
$\\{\mu_{k}\\}$ with $K\rightarrow\infty$. Now let
$\nu_{\ell}^{\prime}=(1-\nu_{\ell})\stackrel{{\scriptstyle
iid}}{{\sim}}\mbox{Beta}(\iota\kappa,\iota)$, so
$\omega_{h}=(1-\nu_{h}^{\prime})\prod_{\ell=1}^{h-1}\nu_{\ell}^{\prime}=\prod_{\ell=1}^{h-1}\nu_{\ell}^{\prime}-\prod_{\ell=1}^{h}\nu_{\ell}^{\prime}$.
By cancelling terms in the cumulative summation
$\pi_{k}=\sum_{h=1}^{k}\omega_{h}$, the result for $1-\pi_{k}$ follows. ∎

######  Lemma 6.

Let $\varepsilon_{n}\to 0$ with
$\varepsilon_{n}^{1/K_{0n}}>\kappa/(\kappa+1)$. For the CUSP prior (5)-(6) and
a positive constant $C>1$, the remainder term $R_{n}=\sum_{k\geq
K_{0n}}\omega_{k}$ satisfies

$\displaystyle\mathbb{P}(R_{n}>\varepsilon_{n})\leq\exp(-CK_{0n}).$ (14)

###### Proof.

Rewriting the reminder as $R_{n}=\prod_{h=1}^{K_{0n}}(1-\nu_{h}),$ the _a
priori_ expectation of the remainder is

$\mathbb{E}(R_{n})=\prod_{h=1}^{K_{0n}}\mathbb{E}(1-\nu_{h})=\left(\frac{\kappa}{\kappa+1}\right)^{K_{0n}}$

from which we conclude that

$\mathbb{P}(R_{0n}>\varepsilon_{n})\leq\frac{1}{\varepsilon_{n}}\left(\frac{\kappa}{\kappa+1}\right)^{K_{0n}}=\left(\frac{1}{\varepsilon_{n}^{1/K_{0n}}}\frac{\kappa}{\kappa+1}\right)^{K_{0n}}\leq\exp(-CK_{0n}),$

where the first inequality follows from Markov’s inequality, and the last one
from the assumptions on the upper bound for $\varepsilon_{n}$. ∎

Before proving Theorem 1, we introduce the following lemma:

######  Lemma 7.

Let $1<r<\sqrt{n}$, $a_{1}>1/2$, and $a_{1}\leq a_{2}$. Under the CUSP and
NMIG in (5)-(6) and (7), respectively, then _a priori_

$\displaystyle\mathbb{P}\left(|\eta_{h}|<\frac{r}{\sqrt{n}}\right)>\left\\{1-\left(\frac{\kappa}{\kappa+1}\right)^{h}\right\\}\left(1-\frac{1}{1+\frac{a_{1}r^{2}}{nv_{0}a_{2}}}\right)$

###### Proof.

First, we consider the conditional probability

$\mathbb{P}\left(|\eta_{h}|<\frac{r}{\sqrt{n}}\mid\pi_{h}\right)$

and then integrate $\pi_{h}$ with respect to the CUSP. We know that, in a
neighborhood of zero, the density function of a Student $t$ random variable
with degrees of freedom greater or equal to 1 is bounded below by the density
of a Cauchy random variable. The cumulative distribution function of a Cauchy
random variable is

$F(t)=\frac{1}{2}+\frac{1}{\pi}\arctan(t).$

Hence, for sufficiently large $n$ we can write

$\displaystyle\mathbb{P}\left(|\eta_{h}|<\frac{r}{\sqrt{n}}\mid\pi_{h}\right)$
$\displaystyle\geq$
$\displaystyle(1-\pi_{h})\frac{2}{\pi}\arctan\left(\frac{r\sqrt{a_{1}}}{\sqrt{a_{2}n}}\right)+\pi_{k}\frac{2}{\pi}\arctan\left(\frac{r\sqrt{a_{1}}}{\sqrt{a_{2}nv_{0}}}\right)$
$\displaystyle\geq$
$\displaystyle\frac{2\pi_{h}}{\pi}\arctan\left(\frac{r\sqrt{a_{1}}}{\sqrt{nv_{0}a_{2}}}\right)$
$\displaystyle\geq$
$\displaystyle\pi_{h}\left(1-\frac{1}{1+\frac{a_{1}r^{2}}{nv_{0}a_{2}}}\right)$

and the last inequality holds since $2/\pi\arctan(t)>(1-1/(1+t^{2}))$ for
$t\in(0,1)$. Finally, integrating out $\pi_{h}$ yields

$\mathbb{P}\left(|\eta_{h}|<\frac{r}{\sqrt{n}}\right)\geq\mathbb{E}\left(\pi_{h}\right)\left(1-\frac{1}{1+\frac{a_{1}r^{2}}{nv_{0}a_{2}}}\right)=\left\\{1-\left(\frac{\kappa}{\kappa+1}\right)^{h}\right\\}\left(1-\frac{1}{1+\frac{a_{1}r^{2}}{na_{2}}}\right).$

∎

######  Theorem 2.

Let $\varepsilon_{n}\to 0$ with $\varepsilon_{n}^{1/K_{0n}}>\kappa/(\kappa+1)$
and assume (C1)-(C3) and $C>2Ae$. For the CUSP (5)-(6) and NMIG (7) priors,
the posterior distribution satisfies

$\lim_{n\to\infty}\mathbb{E}_{0}\left\\{\mathbb{P}\left(R_{n}>\varepsilon_{n}\mid
y_{1},\ldots,y_{n}\right)\right\\}=0.$

###### Proof.

Let $y^{(n)}=(y_{1},\ldots,y_{n})^{\prime}$ denote the observed data. For
$\eta\in\mathbb{R}^{n}$, let $f_{\eta}(\cdot)$ denote the probability density
function of a $N(\eta,I_{n})$ distribution and $f_{\eta_{i}}(\cdot)$ denote
the related univariate marginal. Let $f_{0}(\cdot)$ be the true density, i.e.,
a Gaussian distribution with mean $\eta_{0}$ and identity covariance. Recall
that $\eta_{0i}\neq 0$ for $i\leq K_{0n}$ while the remaining $\eta_{0i}$ for
$i>K_{0n}$ are all null. Let $\Pi(\cdot)$ the CUSP prior measure on the $\eta$
parameters. Let $E_{n}$ be the target event,

$E_{n}=\\{R_{K_{n}}>\varepsilon_{n}\\}.$

Following Castillo and van der Vaart, (2012) and Rockova, (2018), we write

$\mathbb{P}\left(E_{n}\mid{y}^{(n)}\right)=\frac{\int_{E_{n}}\prod_{i=1}^{n}\frac{f_{\eta_{i}}(y_{i})}{f_{0i}(y_{i})}d\Pi(\eta_{i})}{\int\prod_{i=1}^{n}\frac{f_{\eta_{i}}(y_{i})}{f_{0i}(y_{i})}d\Pi(\eta_{i})}=\frac{N_{n}}{D_{n}}.$
(15)

We introduce an event $A_{n}$ with large probability under the true data
generating process. Specifically, let $R=\\{h=K_{n0}+1,\dots,n\\}$ and define

$A_{n}:=\\{D_{n}\geq e^{-r_{n}^{2}}\mathbb{P}(||\eta_{R}||\leq r_{n})\\}.$

If we decompose the probability in (15) in the sum of two complementary
conditional probabilities (conditioning on $A_{n}$ and its complement
$A_{n}^{c}$) we can write

$\mathbb{E}_{0}\left[\mathbb{P}\left(E_{n}\mid{y}^{(n)}\right)\right]\leq\mathbb{E}_{0}\left[\mathbb{P}\left(E_{n}\mid{y}^{(n)}\right)\mathbb{I}_{A_{n}}\right]+\mathbb{P}_{0}(A_{n}^{c})$

Now thanks to Lemma 5.2 of Castillo and van der Vaart, (2012) the second
summand has negligible probability for increasing $n$:
$\mathbb{P}_{0}(A_{n}^{c})\leq\exp\\{-r^{2}_{n}\\}$. Applying (15), we have

$\mathbb{E}_{0}\left[\mathbb{P}\left(E_{n}\mid{y}^{(n)}\right)\right]\leq\frac{\mathbb{P}(E_{n})}{e^{-r_{n}^{2}}\mathbb{P}(||\eta_{R}||\leq
r_{n})}+\exp\\{-r^{2}_{n}\\},$

where $\mathbb{P}(E_{n})$ is simply the prior probability. We can now use
Lemma 6 to bound the numerator and Lemma 7 to bound the denominator.
Specifically, we have $\mathbb{P}(||\eta||\leq
r_{n})\geq\mathbb{P}(|\eta_{K_{0n}}|\leq r_{n}/\sqrt{n})^{n-K_{n0}}$ and
hence, Lemma 7 yields

$\mathbb{P}(||\eta_{R}||\leq
r_{n})\geq\left(1-\left(\frac{\kappa}{\kappa+1}\right)^{K_{0n}}\right)^{n-K_{n0}}\left(1-\frac{1}{1+\frac{a_{1}r_{n}^{2}}{na_{2}v_{0}}}\right)^{n-K_{n0}}.$

Now, we have $\left(\frac{\kappa}{\kappa+1}\right)^{K_{0n}}<AK_{0n}/n$ for
high $n$ and a constant $A>1/2$. Thus, choosing $r^{2}_{n}=K_{0n}$ and
considering (C1), we can write

$\mathbb{P}(||\eta_{R}||\leq
r_{n})\geq\left(1-\frac{AK_{0n}}{n}\right)^{n}\left(1-\frac{1}{n}\right)^{n}\geq\left(\frac{1}{2e}\right)^{AK_{0n}+1},$

where the last inequality follows from $(1-x)^{1/x}>1/(2e)$ for $0<x<0.5$.
Finally, using Lemma 6, we get

$\mathbb{E}_{0}\left[\mathbb{P}\left(E_{n}\mid{y}^{(n)}\right)\right]\leq
2\exp\\{-K_{0n}[C-1-A\log(2e)]+1\\}+\exp\\{-K_{0n}\\}.$

Then, since $C>2Ae>1+A\log(2e)$ for $A>1/2$, we conclude
$\mathbb{E}_{0}\\{\mathbb{P}\left(E_{n}\mid{y}^{(n)}\right)\\}\to 0$. ∎

######  Proposition 7.

For $\varepsilon>0$,
$\mathbb{P}(|\beta_{k,i}|\leq\varepsilon)<\mathbb{P}(|\beta_{k+1,i}|\leq\varepsilon)$
whenever $v_{0}<1$.

###### Proof.

Since $\\{\xi_{k,i}\\}$ are iid (upon marginalization over
$\\{m_{\xi_{k,i}}\\}$), it follows that
$\mathbb{P}(|\beta_{k,i}|\leq\varepsilon)=\mathbb{P}(|\eta_{k}\xi_{k,i}|\leq\varepsilon)=\mathbb{P}(|\eta_{k}\xi_{k+1,i}|\leq\varepsilon)$.
Then we proceed directly:

$\displaystyle\mathbb{P}(|\beta_{k,i}|\leq\varepsilon)$
$\displaystyle=\mathbb{P}(|\eta_{k}\xi_{k+1,i}|\leq\varepsilon)=\mathbb{E}\\{\mathbb{I}(|\eta_{k}||\xi_{k+1,i}|\leq\varepsilon)\\}$
$\displaystyle=\mathbb{E}\big{[}\mathbb{E}\big{\\{}\mathbb{I}(|\eta_{k}||\xi_{k+1,i}|\leq\varepsilon)\big{|}\xi_{k+1,i}\big{\\}}\big{]}$
$\displaystyle=\mathbb{E}\big{[}\mathbb{P}\big{(}|\eta_{k}||\xi_{k+1,i}|\leq\varepsilon\big{|}\xi_{k+1,i}\big{)}\big{]}$
$\displaystyle<\mathbb{E}\big{[}\mathbb{P}\big{(}|\eta_{k+1}||\xi_{k+1,i}|\leq\varepsilon\big{|}\xi_{k+1,i}\big{)}\big{]}$
$\displaystyle=\mathbb{P}(|\beta_{k+1,i}|\leq\varepsilon),$

where the inequality follows from Corollary 4. ∎

######  Proposition 8.

Let $\theta^{(K)}=\\{\theta_{k}\\}_{k=1}^{K}$ denote the sequence
$\\{\theta_{k}\\}_{k=1}^{\infty}$ truncated at $K$. For
$0<v_{0}<\varepsilon<1$, we have
$\mathbb{P}\left\\{d_{\infty}(\theta,\theta^{(K)})>\varepsilon\right\\}\leq\kappa\\{\kappa/(1+\kappa)\\}^{K}$.

###### Proof.

For the prior (7), note that $P_{slab}=\delta_{1}$ and
$P_{spike}=\delta_{v_{0}}$, so
$P_{slab}\\{\bar{\mathbb{B}}_{\varepsilon}(0)\\}=1$ if $\varepsilon<1$ and
$P_{spike}\\{\bar{\mathbb{B}}_{\varepsilon}(0)\\}=0$ if $\varepsilon>v_{0}$,
where $\bar{\mathbb{B}}_{\varepsilon}(0)$ denotes the complement of
$\mathbb{B}_{\varepsilon}(0)$. Using the proof of Proposition 1, we have
$\mathbb{P}\left\\{d_{\infty}(\theta,\theta^{(K)})>\varepsilon\right\\}=\mathbb{P}\left\\{\sup_{k>K}|\theta_{k}|>\varepsilon\right\\}\leq\sum_{k>K}\mathbb{P}\left(|\theta_{k}|>\varepsilon\right)=\sum_{k>K}\\{\kappa/(1+\kappa)\\}^{k}=\kappa\\{\kappa/(1+\kappa)\\}^{K}$.
∎

##  9 MCMC algorithm

We outline the Gibbs sampling algorithm for the semiparametric functional
factor model (SFFM). Specifically, we present full conditional updates for the
model with the following priors on the parametric components:

$[\alpha_{\ell,i}|\sigma_{\alpha_{\ell}}]\stackrel{{\scriptstyle
indep}}{{\sim}}N(0,\sigma_{\alpha_{\ell}}^{2}),\quad\sigma_{\alpha_{\ell}}\stackrel{{\scriptstyle
iid}}{{\sim}}C^{+}(0,1),\quad p(\sigma_{\epsilon}^{2})\propto
1/\sigma_{\epsilon}^{2},$

and assume a generic prior $p(\gamma)$ for the nonlinear parameter $\gamma$
(if unknown). The sampling step for $\gamma$ is omitted if there are no
unknown nonlinear parameters.

### Gibbs sampling algorithm:

  1. 1.

Imputation: for any unobserved $\tau^{*}$ for functional observation $i$,
sample

$[y_{i}(\tau^{*})|-]\stackrel{{\scriptstyle
indep}}{{\sim}}N\Big{(}\sum_{\ell=1}^{L}g_{\ell}(\tau^{*};\gamma)\alpha_{\ell,i}+\sum_{k=1}^{K}f_{k}(\tau^{*})\beta_{k,i},\sigma_{\epsilon}^{2}\Big{)};$

  2. 2.

Nonlinear parameter: if $\gamma$ is unknown, sample

$p(\gamma|-)\propto
p(\gamma)\exp(-\frac{1}{2}\sum_{i=1}^{n}||\bm{y}_{i}-\bm{G}_{\gamma}\bm{\alpha}_{i}-\bm{F}\bm{\beta}_{i}||^{2}/\sigma_{\epsilon}^{2})$

using a slice sampler (Neal,, 2003);

  3. 3.

Nonparametric parameters: for $k=1,\ldots,K$ (in random order),

    1. (a)

Sample the unconstrained coefficients

$\bm{\psi}_{k}^{0}\sim
N\left(\bm{Q}_{\psi_{k}}^{-1}\bm{\ell}_{\psi_{k}},\bm{Q}_{\psi_{k}}^{-1}\right)$

where
$\bm{Q}_{\psi_{k}}=\sigma_{\epsilon}^{-2}(\bm{B}^{\prime}\bm{B})\sum_{i=1}^{n}\beta_{k,i}^{2}+\lambda_{f_{k}}\bm{\Omega}$
and
$\bm{\ell}_{\psi_{k}}=\sigma_{\epsilon}^{-2}\bm{B}^{\prime}\sum_{i=1}^{n}\big{\\{}\beta_{k,i}\big{(}\bm{y}_{i}-\bm{G}_{\gamma}\bm{\alpha}_{i}-\sum_{\ell\neq
k}\bm{f}_{\ell}\beta_{\ell,i}\big{)}\big{\\}}$ (see Kowal and Bourgeois,, 2020
for efficient techniques);

    2. (b)

Compute the constrained vector

$\bm{\psi}_{k}=\bm{\psi}_{k}^{0}-\bm{Q}_{\psi_{k}}^{-1}\bm{C}_{k}^{\prime}\left(\bm{C}_{k}\bm{Q}_{\psi_{k}}^{-1}\bm{C}_{k}^{\prime}\right)^{-1}\bm{C}_{k}\bm{\psi}_{k}^{0}$

where
$\bm{C}_{k}=(\bm{G}_{\gamma},\bm{f}_{1},\ldots,\bm{f}_{k-1},\bm{f}_{k+1},\ldots,\bm{f}_{K})^{\prime}\bm{B}$;

    3. (c)

Normalize: $\beta_{k,i}\rightarrow\beta_{k,i}\|\bm{B}\bm{\psi}_{k}\|$,
$\bm{\psi}_{k}\rightarrow\bm{\psi}_{k}\|\bm{B}\bm{\psi}_{k}\|$, and
$\bm{f}_{k}=\bm{B}\bm{\psi}_{k}$;

    4. (d)

Sample the smoothing parameter

$[\lambda_{f_{k}}|-]\sim\mbox{Gamma}\big{(}(J+1)/2,\bm{\psi}_{k}^{\prime}\bm{\Omega}\bm{\psi}_{k}/2\big{)}\mbox{
truncated to }(10^{-8},\infty)$

where $J$ is the number of columns of $\bm{B}$;

  4. 4.

Ordered spike-and-slab parameters: let
$y_{k,i}^{F}=\bm{f}_{k}^{\prime}\bm{y}_{i}$ for $i=1,\ldots,n$ and
$k=1,\ldots,K$;

    1. (a)

Sample $[m_{\xi_{k,i}}|-]$ from $\\{-1,1\\}$ with
$\mathbb{P}(m_{\xi_{h,i}}=1|-)=1/\\{1+\exp(-2\xi_{h,i})\\}$;

    2. (b)

Sample $[\xi_{k,i}|-]\sim
N(Q_{\xi_{k,i}}^{-1}\ell_{\xi_{k,i}},Q_{\xi_{k,i}}^{-1})$ where
$Q_{\xi_{k,i}}=\eta_{k}^{2}/\sigma_{\epsilon}^{2}+1$ and
$\ell_{\xi_{k,i}}=\eta_{k}y_{k,i}^{F}/\sigma_{\epsilon}^{2}+m_{\xi_{k,i}}$;

    3. (c)

Sample $[\eta_{k}|-]\sim
N(Q_{\eta_{k}}^{-1}\ell_{\eta_{k}},Q_{\eta_{k}}^{-1})$ where
$Q_{\eta_{k}}=\sum_{i=1}^{n}\xi_{k,i}^{2}/\sigma_{\epsilon}^{2}+(\theta_{k}\sigma_{k}^{2})^{-1}$
and
$\ell_{\eta_{k}}=\sum_{i=1}^{n}\xi_{k,i}y_{k,i}^{F}/\sigma_{\epsilon}^{2}$;

    4. (d)

Rescale $\eta_{k}\rightarrow(\sum_{i=1}^{n}|\xi_{k,i}|/n)\eta_{k}$ and
$\bm{\xi}_{k}\rightarrow(n/\sum_{i=1}^{n}|\xi_{k,i}|)\bm{\xi}_{k}$ and update
$\beta_{k,i}=\xi_{k,i}\eta_{k}$;

    5. (e)

Sample
$[\sigma_{k}^{-2}|-]\sim\mbox{Gamma}\left\\{a_{1}+1/2,a_{2}+\eta_{k}^{2}/(2\theta_{k})\right\\}$;

    6. (f)

Sample
$[\nu_{k}|-]\sim\mbox{Beta}(1+\sum_{h=1}^{K}\mathbb{I}\\{z_{h}=k\\},\kappa+\sum_{h=1}^{K}\mathbb{I}\\{z_{h}>k\\})$
for $k=1,\ldots,K-1$ and update
$\omega_{h}=\nu_{h}\prod_{\ell=1}^{h-1}(1-\nu_{\ell})$ and
$\pi_{k}=\sum_{h=1}^{k}\omega_{h}$;

    7. (g)

Sample
$[\kappa|-]\sim\mbox{Gamma}\\{a_{\kappa}+K-1,b_{\kappa}-\sum_{k=1}^{K-1}\log(1-\nu_{k})\\}$;

    8. (h)

Sample $[z_{k}|-]$ from

$\mathbb{P}(z_{k}=h|-)\propto\begin{cases}\omega_{h}t_{2a_{1}}(\eta_{k};0,\sqrt{v_{0}a_{2}/a_{1}})&h\leq
k\\\ \omega_{h}t_{2a_{1}}(\eta_{k};0,\sqrt{a_{2}/a_{1}})&h>k\end{cases}$

where $t_{d}(x;m,s)$ is the density of the $t$-distribution evaluated at $x$
with mean $m$, standard deviation $s$, and degrees of freedom $d$;

    9. (i)

Update $\theta_{k}=1$ if $z_{k}>k$ and $\theta_{k}=v_{0}$ if $z_{k}\leq k$.

  5. 5.

Parametric linear coefficients: sample

$[\bm{\alpha}_{i}|-]\stackrel{{\scriptstyle
indep}}{{\sim}}N\big{(}\bm{Q}_{\alpha_{i}}^{-1}\bm{\ell}_{\alpha_{i}},\bm{Q}_{\alpha_{i}}^{-1}\big{)},\quad
i=1,\ldots,n$

where
$\bm{Q}_{\alpha_{i}}=\sigma_{\epsilon}^{-2}\bm{I}_{L}+\mbox{diag}(\\{\sigma_{\alpha_{\ell}}^{-2}\\}_{\ell=1}^{L})$
and
$\bm{\ell}_{\alpha_{i}}=\sigma_{\epsilon}^{-2}\bm{G}_{\gamma}^{\prime}\bm{y}_{i}$;

  6. 6.

Parametric variance: using a parameter-expansion of the half-Cauchy
distribution (Wand et al.,, 2011),

    1. (a)

Sample
$[\sigma_{\alpha_{\ell}}^{-2}|-,\xi_{\sigma_{\alpha_{\ell}}}]\stackrel{{\scriptstyle
indep}}{{\sim}}\mbox{Gamma}\big{(}(n+1)/2,\sum_{i=1}^{n}\alpha_{\ell,i}^{2}/2+\xi_{\sigma_{\alpha_{\ell}}}\big{)}$
for $\ell=1,\ldots,L$;

    2. (b)

Sample
$[\xi_{\sigma_{\alpha_{\ell}}}|\sigma_{\alpha_{\ell}}]\stackrel{{\scriptstyle
indep}}{{\sim}}\mbox{Gamma}\big{(}1,\sigma_{\alpha_{\ell}}^{-2}+1\big{)}$ for
$\ell=1,\ldots,L$;

  7. 7.

Observation error variance: sample

$[\sigma_{\epsilon}^{-2}|-]\sim\mbox{Gamma}\Big{(}nm/2,\sum_{i=1}^{n}||\bm{y}_{i}-\bm{G}_{\gamma}\bm{\alpha}_{i}-\bm{F}\bm{\beta}_{i}||^{2}/2\Big{)}.$

######  Remark.

Recall that we orthogonalize
$\bm{G}_{\gamma}^{\prime}\bm{G}_{\gamma}=\bm{I}_{L}$ using a QR decomposition
$\bm{G}_{\gamma}^{0}=\bm{Q}_{\gamma}\bm{R}_{\gamma}$ of the initial basis
matrix $\bm{G}_{\gamma}^{0}=(\bm{g}_{1;\gamma},\ldots,\bm{g}_{L;\gamma})$, and
set $\bm{G}_{\gamma}=\bm{Q}_{\gamma}$. When $\gamma$ is unknown and endowed
with a prior distribution, the QR decomposition is incorporated into the
likelihood evaluations of (4) for posterior sampling of $\gamma$. In addition,
the parametric factors may be recovered on the original scale by setting
$\bm{\alpha}_{i}^{0}=\bm{R}_{\gamma}^{-1}\bm{\alpha}_{i}$, which can be
computed draw by draw within the MCMC sampler.

### Hierarchical prior on the nonlinear parameter

For the pinch force data (Section 6.1), we include the hierarchical prior

$\gamma_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}N(\mu_{\gamma},\sigma_{\gamma}^{2}),\quad\mu_{\gamma}\sim
N(0,10),\quad\sigma_{\gamma}\sim C^{+}(0,1)$

which requires minor modifications to the Gibbs sampling algorithm. The
nonlinear parameter is still sampled using the slice sampler (Neal,, 2003),
but with the modified likelihood

$p(\gamma_{i}|-)\propto
p(\gamma_{i}|\mu_{\gamma},\sigma_{\gamma}^{2})\exp(-\frac{1}{2}||\bm{y}_{i}-\bm{G}_{\gamma}\bm{\alpha}_{i}-\bm{F}\bm{\beta}_{i}||^{2}/\sigma_{\epsilon}^{2})$

specific to each $i=1,\ldots,n$, while the (conditional) prior is the
aforementioned Gaussian distribution. The mean parameter is sampled from its
full conditional distribution,

$[\mu_{\gamma}|-]\sim
N\big{(}Q_{\mu_{\gamma}}^{-1}\ell_{\mu_{\gamma}},Q_{\mu_{\gamma}}^{-1}\big{)},$

where $Q_{\mu_{\gamma}}=n\sigma_{\gamma}^{-2}+1/10$ and
$\ell_{\mu_{\gamma}}=\sigma_{\gamma}^{-2}\sum_{i=1}^{n}\gamma_{i}$. The
sampler for the standard deviation parameter uses a half-Cauchy parameter-
expansion (Wand et al.,, 2011):

$\displaystyle[\sigma_{\gamma}^{-2}|-,\xi_{\sigma_{\gamma}}]$
$\displaystyle\sim\mbox{Gamma}\big{(}(n+1)/2,\sum_{i=1}^{n}(\gamma_{i}-\mu_{\gamma})^{2}/2+\xi_{\sigma_{\gamma}}\big{)}$
$\displaystyle[\xi_{\sigma_{\gamma}}|\sigma_{\gamma}]$
$\displaystyle\sim\mbox{Gamma}\big{(}1,\sigma_{\gamma}^{-2}+1\big{)}.$

### Dynamic models for the parametric factors

Since yield curve data (Section 6.2) are time-ordered for $i=1,\ldots,n$, we
specified a dynamic model for the parametric factors

$\alpha_{\ell,i}=\mu_{\ell}+\phi_{\ell}(\alpha_{\ell,i-1}-\mu_{\ell})+\zeta_{\ell,i},\quad\zeta_{\ell,i}\stackrel{{\scriptstyle
indep}}{{\sim}}N(0,\sigma_{\zeta_{\ell}}^{2})$

with the additional priors

$\mu_{\ell}\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,10^{6}),\quad(\phi_{\ell}+1)/2\stackrel{{\scriptstyle
iid}}{{\sim}}\mbox{Beta}(5,2),\quad\sigma_{\zeta_{\ell}}\stackrel{{\scriptstyle
iid}}{{\sim}}C^{+}(0,1)$

for $\ell=1,\ldots,L$. Conditional on $\\{\alpha_{\ell,i}\\}$, the sampling
steps for these parameters are straightforward; see Kowal, (2020) for details.
Note that the prior on the (shifted and scaled) autoregression coefficients
implies that the dynamic factors $\\{\alpha_{\ell,i}\\}$ are stationary, and
therefore so is the functional time series $\\{Y_{i}\\}$.

The dynamic factors can be sampled using state space simulation methods; we
use Durbin and Koopman, (2002) implemented in the KFAS package in R.
Specifically, Corollary 2 implies that the observation equation reduces to

$\bm{g}_{\ell;\gamma}^{\prime}\bm{y}_{i}=\alpha_{\ell,i}+\epsilon_{\ell,i},\quad\epsilon_{\ell,i}\stackrel{{\scriptstyle
indep}}{{\sim}}N(0,\sigma_{\epsilon_{i}}^{2})$

for $\ell=1,\ldots,L$ and $i=1,\ldots,n$, while the evolution equation is
given by the independent AR(1) models for the dynamic factors
$\\{\alpha_{\ell,i}\\}_{i=1}^{n}$ for each $\ell=1,\ldots,L$. Using the state
space construction, we can sample the factors $\\{\alpha_{\ell,i}\\}$ jointly
in $O(nL)$ computational complexity.

The yield curve model also includes a stochastic volatility model on the
observation error variance. The priors and sampling steps are identical to
those in Kowal, (2020).

### Modifications for the Gaussian process alternative

The PFFM+gp (3) was introduced as a seemingly reasonable alternative to (1)
that replaced the nonparametric term with a Gaussian process,
$Y_{i}(\tau)=\sum_{\ell=1}^{L}\alpha_{\ell,i}g_{\ell}(\tau;\gamma)+h_{i}(\tau)$,
where $h_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{GP}(0,\mathcal{K}_{h})$ is a Gaussian process with mean
zero and covariance function $\mathcal{K}_{h}$. Equivalently, each latent
curve $Y_{i}$ is a Gaussian process centered at the parametric term,
$Y_{i}\stackrel{{\scriptstyle
indep}}{{\sim}}\mathcal{GP}\\{\sum_{\ell=1}^{L}\alpha_{\ell,i}g_{\ell}(\cdot;\gamma),\mathcal{K}_{h}\\}$.
Using the observation model (2) with the evaluation points
$\\{\tau_{j}\\}_{j=1}^{m}$, the PFFM+gp modifies (4) to be

$\bm{y}_{i}=\bm{G}_{\gamma}\bm{\alpha}_{i}+\bm{h}_{i}+\bm{\epsilon}_{i},\quad\bm{\epsilon}_{i}\stackrel{{\scriptstyle
indep}}{{\sim}}N(\bm{0},\sigma_{\epsilon_{i}}^{2}\bm{I}_{m})$

where $\bm{h}_{i}=(h_{i}(\tau_{1}),\ldots,h_{i}(\tau_{m}))^{\prime}\sim
N(\bm{0},\bm{K}_{h})$ and
$\bm{K}_{h}=\\{\mathcal{K}_{h}(\tau_{j},\tau_{\ell})\\}_{j,\ell=1}^{m}$.
Informally, each $\bm{h}_{i}$ captures the functional (within-curve)
variability that is unexplained by the parametric term. We implement this
model using a Matérn covariance function with smoothness parameter 2.5, i.e.,
$\mathcal{K}_{h}(\tau_{1},\tau_{2})=\sigma_{h}^{2}[1+\sqrt{5}|\tau_{1}-\tau_{2}|/\rho_{h}+5|\tau_{1}-\tau_{2}|^{2}/(3\rho_{h}^{2})\exp\\{-\sqrt{5}|\tau_{1}-\tau_{2}|/\rho_{h}\\}]$.
We assign the prior $\sigma_{h}\sim C^{+}(0,1)$ and fix $\rho_{h}=0.5$ to
provide smooth paths and reduce the computational burden; the results are not
sensitive to this choice.

To adapt the SFFM Gibbs sampler for the PFFM+gp, we only require modifications
of the sampler for $\\{\bm{\alpha}_{i}\\}$ and sampling steps for the Gaussian
process parameters. Specifically, we sample $\\{\bm{h}_{i}\\}$ after
integrating out the parametric factors $\\{\bm{\alpha}_{i}\\}$, which induces
a joint draw of $\\{\bm{\alpha}_{i},\bm{h}_{i}\\}$ to encourage MCMC
efficiency. Let $\bm{R}_{h}=\sigma_{h}^{-2}\bm{K}_{h}$ denote the (fixed)
correlation matrix. Assuming $\alpha_{\ell,i}\stackrel{{\scriptstyle
indep}}{{\sim}}N(0,\sigma_{\alpha_{\ell}}^{2})$, the requisite (marginalized)
full conditional distributions are $[\bm{h}_{i}\mid-]\stackrel{{\scriptstyle
indep}}{{\sim}}N(\bm{Q}_{h_{i}}^{-1}\bm{\ell}_{h_{i}},\bm{Q}_{h_{i}}^{-1})$
for $i=1,\ldots,n$, where
$\bm{Q}_{h_{i}}=\sigma_{h}^{-2}\bm{R}_{h}^{-1}+\bm{\Sigma}_{h_{i}}^{-1}$,
$\bm{\ell}_{h_{i}}=\bm{\Sigma}_{h_{i}}^{-1}\bm{y}_{i}$, and the marginal
precision
$\bm{\Sigma}_{h_{i}}^{-1}=(\bm{G}_{\gamma}\bm{\Sigma}_{\alpha}\bm{G}_{\gamma}^{\prime}+\sigma_{\epsilon}^{2}\bm{I}_{m})^{-1}=\sigma_{\epsilon}^{-2}\\{\bm{I}_{m}-\bm{G}_{\gamma}(\mbox{diag}\\{\sigma_{\epsilon}^{2}/\sigma_{\alpha_{\ell}}^{2}+1\\}_{\ell=1}^{L})\bm{G}_{\gamma}^{\prime}\\}$
is simplified due to the Woodbury identity and the orthogonality of
$\bm{G}_{\gamma}$. The sampling steps for $\sigma_{h}$ proceed using the same
parameter expansion from Wand et al., (2011) as in Section 9.

We emphasize that _without_ marginalizing over $\\{\bm{\alpha}_{i}\\}$, the
MCMC efficiency for both $\\{\bm{\alpha}_{i}\\}$ and $\\{\bm{h}_{i}\\}$
deteriorates significantly. Hence, our implementation—using both the joint
sampling step and the Woodbury simplifications above—presents a very favorable
sampling algorithm for the PFFM+gp. Nonetheless, both the computational
efficiency and the MCMC efficiency of the PFFM+gp lag far behind those for the
SFFM and the PFFM, both in simulations and real data analyses. Lastly, we note
that the key marginalizations needed for joint sampling of
$\\{\bm{\alpha}_{i},\bm{h}_{i}\\}$ become significantly more challenging when
the model for the parametric factors is more complex (e.g., dynamic); by
comparison, the proposed SFFM avoids this issue entirely.

##  10 Additional simulation results

###  10.1 Linear template: $K_{true}=8$

To assess performance in the presence of many nonparametric factors, we
reproduce the simulation analysis from Section 5 in the case of $K_{true}=8$.
For the proposed SFFM, we increase the bound on the number of factors from
$K=10$ to $K=15$. Figure 8 consolidates the point and interval estimation and
prediction for the linear template; similar results are presented for the
Nelson-Siegel template in Section 10.2. The results are similar to those for
$K_{true}=3$ in the main paper, and confirm the improvements offered by the
SFFM when $K_{true}>0$.

Figure 8: Assessments of point prediction (left) and interval prediction
(right) for $K_{true}=8$. Results are similar to those for $K_{true}=3$.

###  10.2 Nelson-Siegel template

We augment our simulation study to include the Nelson and Siegel, (1987)
template,

$g_{1}(\tau;\gamma)=1,\quad
g_{2}(\tau;\gamma)=\\{1-\exp(-\tau\gamma)\\}/(\tau\gamma),\quad
g_{3}(\tau;\gamma)=g_{2}(\tau;\gamma)-\exp(-\tau\gamma)$

with unknown $\gamma>0$. The synthetic data-generating process proceeds as
before, with following modifications: the true nonlinear parameter is set to
$\gamma^{*}=0.0609$ (Diebold and Li,, 2006; Bianchi et al.,, 2009) and the
nonparametric basis matrix $\bm{F}^{*}$ is orthogonalized to the parametric
basis $\bm{G}_{\gamma^{*}}^{*}$ using a QR-decomposition. As in the yield
curve application, we use a Gamma prior for $\gamma$ with prior mean 0.0609
and prior variance 0.5. The results are given in Figures 9, 10, and 11 and
confirm the results for the linear template in the main paper.

Figure 9: Root mean squared (prediction) errors for $\bm{Y}^{*}$ (top) and
$\\{\alpha_{\ell,i}^{*}\\}$ (bottom) for $K_{true}=0$ (left) and $K_{true}=3$
(right) under the Nelson and Siegel, (1987) template with unknown $\gamma$.
All methods perform similarly for $K_{true}=0$, while the SFFM outperforms all
competitors when $K_{true}>0$. In this setting, the PFFM+gp offers better
prediction than the PFFM, but cannot accurately estimate the parametric
factors.

Figure 10: Mean 95% prediction (top) and credible (bottom) interval widths for
$\bm{y}_{i}^{*}$ and $\\{\alpha_{\ell,i}^{*}\\}$, respectively, with empirical
coverage (annotations) for $K_{true}=0$ (left) and $K_{true}=3$ (right) under
the Nelson and Siegel, (1987) template with unknown $\gamma$. The SFFM
provides narrow interval estimates with (nearly) the correct nominal coverage
in all cases. The PFFM intervals are excessively wide for $K_{true}>0$, while
the PFFM+gp performs poorly for $\\{\alpha_{\ell,i}^{*}\\}$ in all cases.

Figure 11: Probability score $\mathbb{P}(K^{*}=K_{true}\mid\bm{y})$ for
$K_{true}=0$ (left) and $K_{true}=3$ (right) under the Nelson and Siegel,
(1987) template with unknown $\gamma$. The proposed ordered spike-and-slab
prior provides better rank selection and inference than the finite mixture
alternative.

###  10.3 Assessing inclusion probabilities for the nonparametric term

To further evaluate the posterior distribution for the effective number of
nonparametric terms, we report ranked probability scores for
$\mathbb{P}(K^{*}\mid\bm{y})$ under the SFFM and SFFM-fmm in Figure 12. Both
methods perform similarly for $K_{true}=8$, while the proposed SFFM dominates
in the remaining cases. Next, we study the probability of inclusion for _any_
nonparametric terms, $\mathbb{P}(K^{*}>0|\bm{y})$. When $K_{true}>0$, we find
that $\mathbb{P}(K^{*}>0|\bm{y})=1$ in all scenarios, so the SFFMs both offer
excellent power for detecting the presence of a nonparametric addition. For
the more challenging case of the true parametric model with $K_{true}=0$, we
report the summary statistics of $\mathbb{P}(K^{*}>0|\bm{y})$ across 100
simulations in Table 4. The overestimation issues for the finite mixture
alternative are apparent: compared to the SFFM, the SFFM-fmm assigns far
greater inclusion probabilities for the nonparametric term when no such term
is needed.

Figure 12: Negatively-oriented ranked probability score for $\mathbb{P}(K^{*}\mid\bm{y})$ for $K_{true}=0$ (left), $K_{true}=3$ (center), and $K_{true}=8$ (right). These results are for the linear template; results for the Nelson and Siegel, (1987) template with unknown $\gamma$ are nearly identical and are omitted.  Template | Model | Min | 1st quartile | Median | 3rd quartile | Max
---|---|---|---|---|---|---
Linear | SFFM | 0.070 | 0.117 | 0.142 | 0.183 | 0.370
SFFM-fmm | 0.123 | 0.229 | 0.280 | 0.354 | 0.642
Nelson-Siegel | SFFM | 0.075 | 0.118 | 0.138 | 0.180 | 1.000
SFFM-fmm | 0.146 | 0.219 | 0.265 | 0.314 | 1.000
Table 4: Summary statistics of $\mathbb{P}(K^{*}>0|\bm{y})$ across 100
simulations.