# Nonhomothetic CES Preferences:
Closed-Form, Aggregation, and Microfoundation Results

Clement E. Bohr Northwestern University. Email<EMAIL_ADDRESS>Martí Mestieri UPF, CREi, BSE and FRB of Chicago. Email:
<EMAIL_ADDRESS>Emre Enes Yavuz Northwestern University.
<EMAIL_ADDRESS>
# Aggregation and Closed-Form Results
for Nonhomothetic CES Preferences

Clement E. Bohr Northwestern University. Email<EMAIL_ADDRESS>Martí Mestieri UPF, CREi, BSE and FRB of Chicago. Email:
<EMAIL_ADDRESS>Emre Enes Yavuz Northwestern University.
<EMAIL_ADDRESS>
###### Abstract

We provide four novel results for nonhomothetic Constant Elasticity of
Substitution preferences (Hanoch,, 1975). First, we derive a closed-form
representation of the expenditure function of nonhomothetic CES under
relatively flexible distributional assumptions of demand and price
distribution parameters. Second, we characterize aggregate demand from
heterogeneous households in closed-form, assuming that household total
expenditures follow an empirically plausible distribution. Third, we leverage
these results to study the Euler equation arising from standard intertemporal
consumption-saving problems featuring within-period nonhomothetic CES
preferences. Finally, we show that nonhomothetic CES expenditure shares arise
as the solution of a discrete-choice logit problem.

##  1 Introduction

Nonhomothetic Constant Elasticity of Substitution (CES) preferences were
proposed by Hanoch, (1975) and first applied to a general equilibrium setting
by Comin et al., (2021). Nonhomothetic CES preferences have proven useful
since they can be integrated into standard applied general equilibrium models
that rely on the assumption of constant price elasticity of preferences or
production functions. They have been subsequently adopted in applied general
equilibrium exercises in international trade, technology adoption, wage
inequality, business cycles, and to assess the importance of
nonhomotheticities for welfare measurement, e.g., among others, Matsuyama,
(2019), Sposi et al., (2021), Lashkari et al., (2019), Comin et al., (2020),
Jaravel and Lashkari, (2022) and Baqaee and Burstein, (2022).

Nonhomothetic CES preferences (or production functions) are defined through an
implicit equation. This implicit formulation carries over, in general, to the
expenditure function. Moreover, since these preferences do not belong to the
Gorman class, they do not generally admit a representative consumer. This
paper presents a particular case in which these results are overturned. First,
we provide a set of sufficient parametric assumptions such that the
expenditure function admits a closed form. This formulation allows for an
arbitrary correlation between taste parameters and prices faced by a
household. Moreover, it allows for a relatively rich distribution of goods
with different income elasticities—a gamma distribution. The gamma
distribution is flexible enough to generate a density function that is
strictly declining or hump-shaped.111The gamma distribution is generically
skewed, but it also converges to a normal distribution for an arbitrarily
large ”shape parameter,” which we define in the next section.

Second, building on the previous closed-form result, we show that under an
empirically relevant distribution of household total expenditures, we can
characterize the aggregate demand of each good in closed form. To obtain this
result, we assume that households’ total expenditures are distributed
according to the Amoroso distribution. This distribution encompasses many
well-known distributions as particular cases, e.g., log-normal or Pareto, and
it has been shown to fit household income data remarkably well (Salem and
Mount,, 1974). Indeed it was developed and applied by Amoroso himself in 1925
to provide an accurate representation of the income distribution (Amoroso,,
1925).

Third, building on the closed-form representation of the expenditure function,
we explore the Euler equation arising from a standard intertemporal
consumption/savings problem in which the intertemporal utility is CRRA. We
provide a series of results depending on the CRRA coefficient. Finally, we
also provide a derivation of nonhomothetic CES emerging from a discrete-choice
logit model, analogous to that in Anderson et al., (1987) for the homothetic
CES.222See also Trottner, (2022) for an independent derivation in the context
of production functions.

The rest of the note is structured as follows. Section 2 briefly sketches the
general definition of nonhomothetic CES. Section 3 provides the closed-form
derivation of the expenditure function and Section 4 provides the aggregation
results across households with heterogeneous expenditure levels. We provide a
description of the Euler equation under these particular cases in Section 5,
and present the discrete-choice result in Section 6.

##  2 Nonhomothetic CES Preferences

We define utility $U$ over a continuum of goods indexed by $i\in[0,1]$
implicitly through a nonhomothetic CES aggregator

$1=\int_{0}^{1}\left(\Omega_{i}U^{-\varepsilon_{i}}C_{i}\right)^{\frac{\rho-1}{\rho}}di,$
(1)

where $C_{i}>0$ and $\Omega_{i}>0$ denote the consumption and the taste
parameter for good $i$, respectively, $\varepsilon_{i}\geq 0$ governs the
expenditure elasticity for good $i$ and $\rho\in(0,1)\cup(1,\infty)$. In what
follows, we assume $\varepsilon_{i}$ and $\Omega_{i}$ follow a joint
distribution such that Equation (1) yields a $U$ that is unique and quasi-
concave in $\\{C_{i}\\}_{i\in[0,1]}$. Hanoch, (1975) and Comin et al., (2021)
provide sufficient conditions for this to be the case. They also show that the
Hicksian demand for good $i$ implied by maximizing utility $U$ defined in
Equation (1) subject to the budget constraint $\int_{0}^{1}p_{i}C_{i}di=E$ is

$C_{i}=\left(\frac{p_{i}}{E}\right)^{-\rho}\left[\Omega_{i}U^{-\varepsilon_{i}}\right]^{\rho-1},$
(2)

where $p_{i}$ is the price of good $i$ and the total expenditure function is,

$E=\left[\int_{0}^{1}\left(\frac{p_{i}}{\Omega_{i}}U^{\varepsilon_{i}}\right)^{1-\rho}di\right]^{\frac{1}{1-\rho}}.$
(3)

##  3 Closed-Form Solutions

We first show that the mapping between total expenditures and utility given by
the Equation (3) has a closed-form solution under the following two conditions
on the joint distribution of $\varepsilon_{i}$ and $\Omega_{i}$:

###### Assumption 1

The prices and taste parameters have a log-linear relationship with
$\\{\varepsilon_{i}\\}_{i\in[0,1]}$,

$\begin{gathered}\ln p_{i}=\xi_{p}\varepsilon_{i}+\nu_{p},\\\
\ln\Omega_{i}=\xi_{\Omega}\varepsilon_{i}+\nu_{\Omega},\end{gathered}$ (4)

where $\nu_{p}$ and $\nu_{\Omega}$ can follow any distribution as long as the
expectation
$\mathbb{E}\left[\left(\frac{e^{\nu_{p}}}{e^{\nu_{\Omega}}}\right)^{1-\rho}\right]$
exists.

Figure 1: An example of a joint distribution of $\varepsilon_{i}$’s and prices
for which Assumption 1 and 2 hold.

The parameters $\xi_{p}$ and $\xi_{\Omega}$ in Equation (4) govern the
correlations between price, taste parameter, and expenditure elasticity. They
can be an endogenous outcome of a general equilibrium model or exogenously
specified.

###### Assumption 2

$\\{\varepsilon_{i}\\}_{i\in[0,1]}$ are distributed following a gamma
distribution,

$\varepsilon_{i}\sim\text{Gamma}(\alpha,\beta),$ (5)

where $\alpha>0$ and $\beta>0$ are the shape and scale parameters of the gamma
distribution.

The density function is given by,

$f(\varepsilon)=\frac{\varepsilon^{\alpha-1}\exp\left[-\frac{\varepsilon}{\beta}\right]}{\beta^{\alpha}\Gamma(\alpha)}.$
(6)

As a reminder to the reader, when the shape parameter is less than or equal to
one, $\alpha\leq 1$, the distribution has a tail shape, but for $\alpha>1$, it
becomes bell-shaped and a uni-modal distribution. For large values of
$\alpha$, the distribution converges to the normal distribution, while the
case $\alpha=1$ corresponds to the exponential distribution.

###### Result 1

Suppose that Assumptions 1 and 2 hold. Then the mapping between utility and
total expenditure is

$\ln U=\frac{\Upsilon}{1-\rho}-\frac{\Psi}{1-\rho}E^{-\frac{1-\rho}{\alpha}}$
(7)

where,

$\Upsilon=\frac{1}{\beta}-(1-\rho)\left(\xi_{p}-\xi_{\Omega}\right)\in\mathbb{R}\quad,\quad\Psi=\frac{\mathcal{M}^{\frac{1}{\alpha}}}{\beta}\in\mathbb{R}_{+}\quad,\quad\mathcal{M}\equiv\mathbb{E}\left[\left(\frac{e^{\nu_{p}}}{e^{\nu_{\Omega}}}\right)^{1-\rho}\right].$

Moreover, consumption choices $C_{i}$ in equation (2) are invariant to the
scaling of $\beta$. Thus, we can normalize $\beta=1$ without loss of
generality.

The closed-form result in Equation (7) follows from direct integration of the
expenditure function.333 After specifying the distribution of $p_{i}$ and
$\Omega_{i}$, Equation 3 becomes  $\displaystyle E^{1-\rho}$
$\displaystyle=\int_{\nu_{\Omega}}\int_{\nu_{p}}\int_{0}^{\infty}\left(\frac{e^{\nu_{p}}}{e^{\nu_{\Omega}}}\right)^{1-\rho}\exp\left[(\xi_{p}-\xi_{\Omega}+\ln
U)\varepsilon)\right]^{1-\rho}f(\varepsilon)f(\nu_{p})f(\nu_{\Omega})=\mathbb{E}\left[\left(\frac{e^{\nu_{p}}}{e^{\nu_{\Omega}}}\right)^{1-\rho}\right]\int_{0}^{\infty}\exp\left[(\xi_{p}-\xi_{\Omega}+\ln
U)\varepsilon)\right]^{1-\rho}f(\varepsilon)$
$\displaystyle=\mathbb{E}\left[\left(\frac{e^{\nu_{p}}}{e^{\nu_{\Omega}}}\right)^{1-\rho}\right]\left(\frac{\tilde{\beta}}{\beta}\right)^{\alpha}\int_{0}^{\infty}\frac{\varepsilon^{\alpha-1}\exp\left[-\frac{\varepsilon}{\tilde{\beta}}\right]}{\tilde{\beta}^{\alpha}\Gamma(\alpha)}=\mathbb{E}\left[\left(\frac{e^{\nu_{p}}}{e^{\nu_{\Omega}}}\right)^{1-\rho}\right]\left(\frac{\tilde{\beta}}{\beta}\right)^{\alpha},$
(8) where
$\tilde{\beta}=\left(\frac{1}{\beta}-(1-\rho)(\xi_{p}-\xi_{\Omega}+\ln
U)\right)^{-1}$ and the last step follows from the integral over a Gamma
distribution density with parameters of $\alpha$ and $\tilde{\beta}$. This
step requires $\tilde{\beta}>0$, and it is possible always to ensure it. Using
the closed-form mapping between the utility level and total expenditure in
Equation (7) and the demand Equation (2), we find that the expenditure
elasticity of good $i$ is

$\eta_{i}=\frac{\partial\ln C_{i}}{\partial\ln
E}=\rho+(1-\rho)\mathcal{M}^{\frac{\beta}{\bar{\varepsilon}}}E^{\beta\frac{\sigma-1}{\bar{\varepsilon}}}\frac{\varepsilon_{i}}{\bar{\varepsilon}}.$
(9)

where $\bar{\varepsilon}=\int_{0}^{1}s_{i}\varepsilon_{i}di$ is the
expenditure-share weighted mean of $\varepsilon_{i}$ and $s_{i}$ denotes the
expenditure share in good $i$.

Next, we show that Equation (9) is homogeneous of degree zero in $\beta$. This
implies that expenditure elasticities and consumption choices are invariant to
the choice of $\beta$. Suppose we scale $\beta$ proportionally by a factor
$k>0$. Using the definition of utility, we find that this is equivalent to
multiplying all $\varepsilon_{i}$’s in Equation (1) with the same constant.
More precisely, we have that if $\varepsilon_{i}$ is distributed gamma (i.e.,
Equation (5) in condition 2 holds) then

$k\varepsilon_{i}\sim\text{Gamma}(\alpha,k\beta).$ (10)

Since the expenditure elasticity of good $i$ is pinned down by the value of
$\varepsilon_{i}$ relative to its mean value
$\overline{\varepsilon}_{i}\equiv\mathbb{E}[\varepsilon_{i}]=\alpha\beta$, the
ratio of the two is independent of $k$. Therefore, $\beta$ can be normalized
without losing generality.

We conclude this section by noting that when the goods are complements, that
is, $\rho<1$, the goods with higher values of $\varepsilon_{i}$ are more
expenditure elastic, and the expenditure elasticity (9) is always positive for
all goods at any expenditure level. When the goods are substitutes, the
relationship is reversed; higher $\varepsilon_{i}$ means less expenditure
elastic goods.444Moreover, the expenditure elasticity becomes negative for a
good $i$ when the expenditure level exceeds the expenditure level
$E^{*}_{i}=\left(\frac{\rho-1}{\rho}\mathcal{M}^{\frac{\beta}{\bar{\varepsilon}}}\frac{\varepsilon_{i}}{\bar{\varepsilon}}\right)^{-\beta\frac{\bar{\varepsilon}}{(\sigma-1)}}$.
As the expenditure level increases, consumption increases relatively more on
goods with a high expenditure elasticity (low $\varepsilon_{i}$). Since the
gamma distribution is defined on the positive real line, $\varepsilon\geq 0$,
eventually, the expenditure share of goods with $\varepsilon>0$ has to
decline.

##  4 Demand Aggregation of Households with Heterogeneous Total Expenditure
Levels

We consider an economy populated by households with heterogeneous expenditures
satisfying Result 1, so that they have an expenditure function given by
Equation (7). In this case, we show that if total expenditure levels follow an
Amoroso distribution, the aggregate demand has a closed-form solution. The
Amoroso distribution is named after the Italian economist Luigi Amoroso who
used it to study income distributions Amoroso, (1925). A central motivation
for Amoroso was to address the counterfactual feature of the Pareto
distribution of having a mode at the lower bound of its support. This
distribution appears to provide a reasonably good representation of total
household expenditures and income (Salem and Mount,, 1974).

Consider household $h$ with total expenditure $E_{h}$ given by Equation (7).
In this case, the expenditure share in good $i$ is

$s_{ih}=\exp(\varepsilon_{i}\Upsilon)\left(\frac{\Omega_{i}}{p_{i}}\right)^{\rho-1}E_{h}^{\rho-1}\exp\left(-\varepsilon_{i}\Psi
E_{h}^{\frac{\rho-1}{\alpha}}\right).$ (11)

We assume the expenditure distribution follows an Amoroso distribution, whose
probability density function is given by,

$f_{E_{h}}(x\mid
l,k,m,n)=\frac{1}{\Gamma(m)}\left|\frac{n}{k}\right|\left(\frac{x-l}{k}\right)^{mn-1}\exp\left\\{-\left(\frac{x-l}{k}\right)^{n}\right\\}.$
(12)

with $m>0$ and $k,l,n\in\mathbbm{R}$. We shall focus on the economically
relevant case $k>0$. In this case, the support of the distribution is $x\geq
l$. The Amoroso distribution has two shape parameters, $m$ and $n$, a scale
parameter $k$, and a location parameter $l$. Since the units of measurement of
expenditures are a free degree of freedom, we focus on the particular
normalization of $l=0$ without loss of generality. It is readily verified by
direct substitution that the Amoroso distribution in Equation (12) has as
particular cases many familiar distributions such as the exponential, gamma,
Fréchet, Weibull distributions, and as limiting cases as $n\rightarrow 0$, the
lognormal and Pareto distributions (depending on what we assume about $m$).555
Over fifty distinct distributions occur as special cases or limiting forms of
the Amoroso distribution; see Crooks, (2015).

###### Assumption 3

The second shape parameter $n$ of the Amoroso distribution (12) satisfies
$n=\frac{\rho-1}{\alpha}$. The location parameter is normalized to zero, $l=0$
without loss of generality.

Figure 2: Amoroso distribution with different shape parameters.

This assumption links the second shape parameter of the Amoroso distribution
to the elasticity of substitution parameter $\rho$ and one of the parameters
of the Gamma distribution of $\varepsilon_{i}$ appearing in Assumption 2. So,
when this assumption holds, we have only two degrees of freedom left in the
Amoroso distribution, $k$ and $m$. Fortunately, varying these two parameters
yields a wide range of potential shapes that can be bell-shaped or
monotonically decreasing (e.g., Weibull, Fréchet or Pareto). Figure 2 also
shows that the changes in the distribution induced by the normalized parameter
$n$ are qualitatively similar to those of $m$ and $k$.

###### Result 2

Consider an economy populated by a continuum of agents indexed by $h$, whose
total expenditure, $E_{h}$, is distributed Amoroso (Equation 12). If
Assumptions 1, 2 and 3 hold, total expenditure shares of good $i$ in the
population are

$s_{i}=\int
s_{ih}f_{E_{h}}(E_{h})dE_{h}=\exp(\varepsilon_{i}\Upsilon)\left(\frac{\Omega_{i}}{p_{i}}\right)^{\rho-1}\frac{\Gamma(m+\alpha)}{\Gamma(m)}\frac{k^{\rho-1}}{\left[1+\varepsilon_{i}\Psi
k^{\frac{\rho-1}{\alpha}}\right]^{m+\alpha}}$ (13)

The expenditure share in terms of the average household expenditure,
$\bar{E}_{h}=k\frac{\Gamma\left(m+\frac{\alpha}{\rho-1}\right)}{\Gamma(m)}$,
is

$s_{i}=\int
s_{ih}f_{E_{h}}(E_{h})dE_{h}=\exp(\varepsilon_{i}\Upsilon)\left(\frac{\Omega_{i}}{p_{i}}\right)^{\rho-1}\frac{\Gamma(m+\alpha)}{\Gamma(m+\frac{\alpha}{\rho-1})}\frac{k^{\rho-2}}{\left[1+\varepsilon_{i}\Psi
k^{\frac{\rho-1}{\alpha}}\right]^{m+\alpha}}\bar{E}_{h}$ (14)

which approximately corresponds to666The result follows from the approximation
of the ratio of gamma functions $\Gamma(m+\alpha/(\rho-1))/\Gamma(m)\approx
m^{\alpha}/(\rho-1)$, which implies
$\Gamma(m+\alpha/(\rho-1))/\Gamma(m)\approx\left(\Gamma(m+\alpha)/\Gamma(m)\right)^{1/(\rho-1)}$.
This approximation is better for $m\gg\alpha/(\rho-1)$, see Frame, (1949).

$s_{i}\approx\exp(\varepsilon_{i}\Upsilon)\left(\frac{\Omega_{i}}{p_{i}}\right)^{\rho-1}\frac{\bar{E}_{h}^{\rho-1}}{\left[1+\varepsilon_{i}\Psi
k^{\frac{\rho-1}{\alpha}}\right]^{m+\alpha}}.$ (15)

Equations (14) and (15) show that the distribution of expenditures affects
aggregate sectoral consumption—not surprisingly, since nonhomothetic
preferences are not Gorman-aggregable. Perhaps interestingly, Equation (14)
implies that aggregate consumption is log separable in the consumption choice
of the household with average expenditures and another term that captures the
importance of inequality in the form of higher moments of the expenditure
distribution.

##  5 Intertemporal Substitution and Euler Equation

Nonhomothetic CES preferences also affect the consumption/savings decision of
households. The intertemporal substitution is affected because households
account for the change in their consumption basket due to changes in their
utility level when they consider how to substitute expenditure between time
periods. In this note, we focus on the case in which the per-period flow
utility is given by a CRRA function,

$v(U)=\frac{U^{1-\theta}-1}{1-\theta}$ (16)

where $\theta$ controls the intertemporal elasticity of substitution. The
household maximizes the discounted sum of flow utility,

$\max_{\\{E_{t},A_{t}\\}_{t=0}^{\infty}}\sum_{t=0}^{\infty}\beta^{t}v(U)$ (17)

subject to within-period preferences being the nonhomothetic CES preferences
defined in Section 2 and the intertemporal budget constraint

$E_{t}+A_{t+1}\leq(1+r_{t})A_{t}+Y_{t},$ (18)

where $A_{t}$ is a risk-free asset in zero net supply ($A_{0}=0$) and $Y_{t}$
is household income.

###### Result 3

  1. (i)

The Euler equation is given by,

$\left(\frac{E_{t+1}}{E_{t}}\right)^{\theta}=\beta(1+r_{t})\left(\frac{\overline{\varepsilon}_{t+1}}{\overline{\varepsilon}_{t}}\right)^{-1}\left(\frac{P_{t+1}}{P_{t}}\right)^{\theta-1}$
(19)

where $P_{t}$ is the ideal price index defined as $\frac{E_{t}}{U_{t}}$ and
$\overline{\varepsilon}_{t}$ is a weighted average $\varepsilon_{i}$ at time
$t$, given by $\overline{\varepsilon}_{t}=\int_{0}^{1}s_{i}\varepsilon_{i}di$
with $s_{i}=\frac{p_{it}C_{it}}{E_{t}}$ being expenditure share of good $i$.

  2. (ii)

If $p_{i}$ and $\Omega_{i}$ follow distributions with Assumption 1 holding,
and hence with a closed form mapping between $U_{t}$ and $E_{t}$ given by
Equation (7), then we have
$\overline{\varepsilon}_{t}=\frac{\alpha}{\Psi}E_{t}^{\frac{1-\rho}{\alpha}}$.
In addition, if we normalize the ideal price index to one, $P_{t}=1$, for all
$t$ and assume $\theta+\frac{1-\rho}{\alpha}\geq 0$, then the Euler equation
becomes,777 Without normalizing the ideal price index, the Euler equation is
given by
$\left(\frac{E_{t+1}}{E_{t}}\right)^{1+\frac{1-\rho}{\alpha}}=\beta(1+r_{t})\exp\left[\frac{(\theta-1)\Psi}{1-\rho}\left(E_{t}^{-\frac{1-\rho}{\alpha}}-E_{t+1}^{-\frac{1-\rho}{\alpha}}\right)\right].$
(20)

$\left(\frac{E_{t+1}}{E_{t}}\right)^{\theta+\frac{1-\rho}{\alpha}}=\beta(1+r_{t}).$
(21)

  3. (iii)

If there are households with identical preferences and their total
expenditures distributed Amoroso at time $t$, $E_{ht}\sim Amoroso(l,k,m,n)$
and the previous result holds, then the expenditure distribution at time $t+1$
is also Amoroso, $E_{ht+1}\sim Amoroso(l,A\times k,m,n)$ where
$A=\left[\beta(1+r_{t+1})\right]^{\frac{\alpha}{\alpha\theta+1-\rho}}$.

The solution of the intertemporal optimization defined in the Equation (17) is
given by

$\frac{v^{\prime}(U_{t})\frac{\partial U_{t}}{\partial
E_{t}}}{v^{\prime}(U_{t+1})\frac{\partial U_{t+1}}{\partial
E_{t+1}}}=\beta(1+r_{t+1}).$ (22)

The household arranges its expenditure between periods so that the ratio of
flow utility $v(U_{t})$ that one additional dollar can buy equalizes to the
discounted rate of return. The curvature driven by CRRA preferences is
standard $v^{\prime}(U_{t})=U_{t}^{-\theta}$. However, $\frac{\partial
U_{t}}{\partial E_{t}}$ is not constant due to preferences being
nonhomothetic. In particular,

$\frac{\partial U_{t}}{\partial
E_{t}}=\frac{\overline{\varepsilon}_{t}^{-1}}{P_{t}}.$ (23)

By Walras’ law, we can normalize $P_{t}=1$. Then, the derivative in (23) is
only determined by the average $\varepsilon_{i}$ weighted by the expenditure
shares. If the goods are complements, $\rho<1$, then the higher
$\varepsilon_{i}$ means more expenditure elastic goods, so
$\overline{\varepsilon}_{t}$ can be thought as an average luxury index of the
consumption basket at a given time. Therefore, in this case, as the household
gets richer and spends relatively more on more luxurious goods, one more
dollar buys less utility. If Assumption 1 holds,
$\overline{\varepsilon}_{t}\propto E_{t}^{\frac{1-\rho}{\alpha}}$. Therefore
the nonhomotheticity acts as if we had some additional curvature in
intertemporal preferences (or lesser curvature if the goods are substitutes,
$\rho>1$).

##  6 Logit Microfoundation of Preferences

In this last section, we show that it is possible to derive the expenditure
shares of the nonhomothetic CES as the probability choice of a logit model
analogously to Anderson et al., (1987). This provides an alternative
aggregation result for nonhomothetic CES. Assume an economy populated by a
mass of households with linear preferences over goods indexed by $i\in[0,1]$.
Households optimally consume a single good according to their indirect utility
function,

$V_{i}=\ln\Omega_{i}+(1-\varepsilon_{i})\ln\frac{E}{p}-\ln\frac{p_{i}}{p}+\mu\nu_{i}$
(24)

where $\frac{E}{p}=U$ denotes the level of "real expenditures" according to
the optimal utility-based price-index $p$, $\nu_{i}$ is an unobservable
preference heterogeneity shock following a standard Gumbel
distribution,888That is, with cumulative distribution is $\exp(\exp(-\nu))$.
and $\mu$ is a parameter governing the degree of preference heterogeneity.
Lastly, $\Omega_{i}$ captures observable preference heterogeneity or tastes.
It is well known that the probability of a household choosing to consume good
$i$ takes a logit form. Under this utility specification, it becomes

$\operatorname{Pr}(i)=\frac{\left(\frac{E}{p}\right)^{\frac{1-\varepsilon_{i}}{\mu}}\left(\Omega_{i}^{-1}\frac{p_{i}}{p}\right)^{-\frac{1}{\mu}}}{\int\left(\frac{E}{p}\right)^{\frac{1-\varepsilon_{i}}{\mu}}\left(\Omega_{i}^{-1}\frac{p_{i}}{p}\right)^{-\frac{1}{\mu}}di}.$
(25)

We set $\mu=-\frac{1}{1-\rho}$ so that the denominator equals one by the
definition of the utility. Furthermore, if all households have the same
expenditure level, $E$, then the demand share of good $i$ is given by the
probability above,

$s_{i}=\frac{p_{i}C_{i}}{E}=\left(\Omega_{\varepsilon
i}^{-1}\frac{p_{i}}{p}\right)^{1-\rho}\left(\frac{E}{p}\right)^{(\varepsilon_{i}-1)(1-\rho)},$
(26)

which is the same as the demand share implied by the nonhomothetic CES
aggregator if we choose
$p=\left[\int\left(\frac{p_{i}}{\Omega_{i}}U^{\varepsilon_{i}-1}\right)^{1-\rho}di\right]^{\frac{1}{1-\rho}}$.
We have shown thus the following result.

###### Result 4

Suppose that indirect utility is given by Equation (24) with
$E=\left[\int\left(\frac{p_{i}}{\Omega_{i}}U^{\varepsilon_{i}}\right)^{1-\rho}di\right]^{\frac{1}{1-\rho}}$.
Consider a population of mass 1 and expenditure $E$ in which individual $i$
maximizes $V_{i}$ given an independent draw of $v_{i}$. Aggregate expenditure
shares in the economy coincide with those of a consumer with nonhomothetic CES
yielding expenditure shares as in Equation (26).

##  7 Conclusion

This note has presented sufficient conditions for a closed-form representation
of the expenditure function of nonhomothetic CES and aggregation of demand of
households with heterogeneous levels of total expenditure. We also analyzed
the resulting Euler equation under this parameterization of nonhomothetic CES
and provided a link with a discrete-choice model. We believe that these
results may be of use in applications that rely on nonhomothetic CES
formulation of preferences or production. Finally, we note that an analogous
exercise can be performed to heterothetic Cobb-Douglas preferences (Bohr et
al.,, 2023).

## References

  * Amoroso, (1925) Amoroso, L. (1925).  Ricerche intorno alla curva dei redditi.  Annali di Matematica Pura ed Applicata, 2:123–159.
  * Anderson et al., (1987) Anderson, S. P., De Palma, A., and Thisse, J.-F. (1987).  The ces is a discrete choice model?  Economics Letters, 24(2):139–140.
  * Baqaee and Burstein, (2022) Baqaee, D. R. and Burstein, A. (2022).  Welfare and Output With Income Effects and Taste Shocks*.  The Quarterly Journal of Economics, 138(2):769–834.
  * Bohr et al., (2023) Bohr, C. E., Mestieri, M., and Robert-Nicoud, F. (2023).  Heterothetic cobb douglas: Theory and applications.  Technical report.
  * Comin et al., (2020) Comin, D., Danieli, A., and Mestieri, M. (2020).  Income-Driven Labor-Market Polarization.  Working Papers 2020-050, Human Capital and Economic Opportunity Working Group.
  * Comin et al., (2021) Comin, D., Lashkari, D., and Mestieri, M. (2021).  Structural Change With Long-Run Income and Price Effects.  Econometrica, 89(1):311–374.
  * Crooks, (2015) Crooks, G. E. (2015).  The amoroso distribution.
  * Frame, (1949) Frame, J. S. (1949).  An approximation to the quotient of gamma function.  The American Mathematical Monthly, 56(8):529–535.
  * Hanoch, (1975) Hanoch, G. (1975).  Production and Demand Models with Direct or Indirect Implicit Additivity.  Econometrica, 43(3):395–419.
  * Jaravel and Lashkari, (2022) Jaravel, X. and Lashkari, D. (2022).  Nonparametric measurement of long-run growth in consumer welfare.  CEP Discussion Papers dp1859, Centre for Economic Performance, LSE.
  * Lashkari et al., (2019) Lashkari, D., Bauer, A., and Boussard, J. (2019).  Information Technology and Returns to Scale.  2019 Meeting Papers 1380, Society for Economic Dynamics.
  * Matsuyama, (2019) Matsuyama, K. (2019).  Engel’s Law in the Global Economy: Demand-Induced Patterns of Structural Change, Innovation, and Trade.  Econometrica, 87(2):497–528.
  * Salem and Mount, (1974) Salem, A. B. Z. and Mount, T. D. (1974).  A convenient descriptive model of income distribution: The gamma density.  Econometrica, 42(6):1115–27.
  * Sposi et al., (2021) Sposi, M., Yi, K.-M., and Zhang, J. (2021).  Deindustrialization and Industry Polarization.  NBER Working Papers 29483, National Bureau of Economic Research, Inc.
  * Trottner, (2022) Trottner, F. (2022).  Who gains from scale?  Technical report.


# Generalizability Under Sensor Failure:
Tokenization + Transformers Enable More Robust Latent Spaces

Geeling Chau, Yujin An11footnotemark: 1, Ahamed Raffey Iqbal11footnotemark: 1,
Soon-Jo Chung, Yisong Yue, Sabera Talukder
California Institute of Technology, Pasadena, CA
<EMAIL_ADDRESS>Equal Contribution

###### Abstract

A major goal in neuroscience is to discover neural data representations that
generalize. This goal is challenged by variability along recording sessions
(e.g. environment), subjects (e.g. varying neural structures), and sensors
(e.g. sensor noise), among others. Recent work has begun to address
generalization across sessions and subjects, but few study robustness to
sensor failure which is highly prevalent in neuroscience experiments. In order
to address these generalizability dimensions we first collect our own
electroencephalography dataset with numerous sessions, subjects, and sensors,
then study two time series models: EEGNet (Lawhern et al., 2018) and TOTEM
(Talukder et al., 2024). EEGNet is a widely used convolutional neural network,
while TOTEM is a discrete time series tokenizer and transformer model. We find
that TOTEM outperforms or matches EEGNet across all generalizability cases.
Finally through analysis of TOTEM’s latent codebook we observe that
tokenization enables generalization.

##  1 Introduction

Neuroscience experiments vary across numerous dimensions including sessions,
subjects and senors (Gonschorek et al., 2021; Saha & Baumert, 2020; Parvizi &
Kastner, 2018). Given this inherent variability, models with strong
generalization properties are desirable. Model generalizability refers to a
model’s zero shot capabilities, or the model’s ability to operate on test
datasets unseen at training time. Prior work studies generalizability along
the datasets’ session and subject dimensions (Peterson et al., 2021; Talukder
et al., 2022; Krumpe et al., 2017), but few study model generalizability under
sensor variability. Sensors vary primarily because of sensor failure and
sensor count differences across recording sessions. Common practice is to
reduce the train and test sets to the intersection of clean and available
sensors then train models. However, this throws away difficult-to-obtain
neural data. It is therefore valuable to build models that can train on all
available data and infer under any type of sensor variability.

To probe generalizability across sessions, subjects, and sensors, we
systematically study two time series models: EEGNet (Lawhern et al., 2018) and
TOTEM (Talukder et al., 2024). EEGNet is a popular convolutional neural
network upon which many other models are built or compared against (Peterson
et al., 2021; Xu et al., 2021). EEGNet intakes a Sensor x Time array and
applies convolutional kernels along the temporal and spatial dimensions, a
common approach inspired by Filter-Bank Common Spatial Patterns (FBCSP) Ang et
al. (2012). TOTEM first learns a sensor-agnostic set of discrete tokens via a
self-supervised vector quantized variational autoencoder, then uses these
pretrained tokens as the input to a transformer classifier. TOTEM exhibits
strong generalization on numerous time series datasets (Talukder et al.,
2024). To study these models’ generalizability across experimental conditions
we create a taxonomy of generalization cases that encompass session, subject,
and sensor variability. We collect a rich electroencephalography dataset that
permits testing of these generalization cases. See Figure 1 for a
visualization of cross subject generalizability under sensor failure. Finally
we explore TOTEM’s latent discrete codebooks, as TOTEM demonstrates the best
performance across experiments in our generalization taxonomy.

Figure 1: Overview. TOTEM and EEGNet train on data from subject A1, which has
no failed sensors. Both TOTEM and EEGNet are then tested on subject B2 with
artificially failed sensors. This is an example of cross subject
generalizability under sensor failure.

##  2 Metholodology

###  2.1 Experimental design

We first define a taxonomy of a baseline and generalization cases, the cases
are as follows: Baseline Case: _Within Session_ \- train and test on the same
session. Generalization Cases: (1) _Cross Session_ \- train on one session,
and test on a separate session from the same subject. (2) _Cross Subject_ \-
train on one subject, and test on a separate subject. (3) _Sensor Failure_ \-
randomly fail X% of sensors, $X\in\\{0,10,20,...,100\\}$. Notably the within
session, cross session, and cross subject cases can be combined with the
sensor failure condition, e.g. one can study cross subject generalizability
under sensor failure (Figure 1). For a visualization of the baseline and
generalizability cases see Figure 2(a).

We then design and collect an electroencepholography (EEG) dataset featuring
high sensor count (128 electrodes) and high trial count (600 trials/session)
across two human subjects. In total this leads to four sub-datasets across
subjects A,B and sessions 1,2: A1, A2, B1, B2. Our 128 sensor setup introduces
relatively large datasets as previous publicly available EEG datasets rarely
exceed 64 sensors (Stieger et al., 2021; Kaya et al., 2018; Tangermann et al.,
2012). Additionally, much prior work with more than 64 channels collects <200
trials per session (Gwon et al., 2023; Xu et al., 2020) we collect 150 trials
per direction (600 trials per session). Designing our experiments to generate
relatively large datasets enables our study of generalization across numerous
experimental conditions.

For each session a subject sat in front of a monitor and fixated on a center
point that randomly changed to $\blacktriangleleft$, $\blacktriangleright$,
$\blacktriangle$, and $\blacktriangledown$; $\blacktriangleleft$ is left hand
movement (HM), $\blacktriangleright$ is right HM, $\blacktriangle$ is both HM,
and $\blacktriangledown$ is both feet movement. In each session there were 150
trials for each direction, lasting 3 seconds each. After data collection we
performed minimal data processing by (1) downsampling from 4096Hz to 2048Hz,
(2) high-pass filtering at 0.5Hz, (3) average referencing, and (4)
standardizing across the entire recording.

We chose to use EEG due to its high variance along these generalization cases,
allowing us to study generalizability in one of the most difficult data
modalities. However, challenges regarding sensor variation are highly
prominent in larger signal-to-noise ratio recording modalities such as
electrocorticography or multi-unit probes (Neuropixels, Utah arrays, DBS
probes). In the future we aim to replicate our EEG experiments with other
neural data recording modalities, see Section 4 for further discussion.

(a)

(b)

Figure 2: (a) Visualization of baseline (within session) and generalizability
cases (cross subject, cross session, sensor failure) illustrated. (b)
Experimental setup. We test on 2 human subjects each with 128 electrodes
generating 600 trials per session. Adapted from García-Murillo et al. (2023)

###  2.2 Models

EEGNet from Lawhern et al. (2018) is a popularly used convolutional neural
network (CNN) for EEG decoding. It consists of temporal and spatial CNNs that
learn kernels swept across these dimensions. Its design mimics a performant
EEG processing pipeline, FBCSP (Ang et al. (2012)), learning temporal kernels
that extract frequency features which then get weighted by learned spatial
filters (Figure 3). EEGNet takes a fixed [Sensors x Time] matrix as input.

Figure 3: EEGNet architecture. Adapted from Lawhern et al. (2018)

TOTEM from Talukder et al. (2024) consists of a sensor-agnostic tokenizer
which uses a learned latent codebook, followed by a time-length flexible and
sensor-count flexible transformer classifier (Figure 4). TOTEM first learns
the latent codebook via self-supervision on the time series signals (Figure
4(a)). Its goal is to learn a set of small time tokens that help reconstruct
the time series signals with lowest MSE. This can be seen as learning a
discrete temporal “vocabulary” that can be used to express the data modality’s
distribution of temporal activity. This tokeniation technique is inspired by
the VQ-VAE (Aaron et al. (2017)) which has been studied to be beneficial in
helping (1) reduce noise in sequential representations and (2) allow
transformers to operate on a meaningful and defined vocabulary set (Talukder
et al., 2024). After the codebook training has converged, it is then frozen
and used to translate the time series inputs into tokenized inputs for the
downstream temporal and spatial transformers (Figure 4(b)). These transformers
are able to operate on different sequence lengths (Vaswani et al. (2017)),
providing benefits for modeling trials with different lengths and sensor
availabilities. We here only study the case where the input sizes are not
changed in order to be able to compare with EEGNet.

Figure 4: TOTEM architecture and training. (a) Learn latent codebok via self-
supervision. (b) Train transformers on tokenized data created from frozen
codebook. Adapted from Talukder et al. (2024)

###  2.3 Training

The 600 trials in each session are broken down into fixed training (80%),
validation (10%), and test (10%) sets, within each of which are then broken
down into 250ms long individual mini-trials. We train and test on these mini-
trials. In the baseline within session case, we train, val, and test on the
same session. In generalization cross session and subject cases, we train and
validate on the first session, then test on the second session. In sensor
failure cases, we train and validate on the original sensors and test on a set
with artificially failed sensors.

Sensor Failure: We simulate failure by randomly zeroing-out X% of test set
sensors, where X$\in\\{0,10,20,...,100\\}$ (11 cases). We cannot remove
sensors as EEGNet requires fixed input size, so we zero-out sensors for both
models. For fair comparison, we ensure both models have the same failed
sensors for each random seed. We study the effect on classification
performance across within session and generalizability cases.

Hyperparameter selection: For each model, we selected a set of hyperparameters
that allowed the models to converge and perform optimally on within session
performance. We fixed these hyperparameters and ran all our modeling
experiments with the same parameters across datasets.

##  3 Results & Discussion

We report decoding accuracy across all generalization cases (Figure 5) and
analyze TOTEM’s latent spaces to understand how its representations generalize
(Figure 6).

###  3.1 Decoding Accuracy Comparison

To evaluate the robustness of the latent spaces learned by each of the models,
we tested their performance on the baseline case (within session) and zero-
shot to generalization cases (cross session, cross subject, and across sensor
failure levels).

Figure 5: Classifier performance under all generalizability cases. (a) TOTEM
vs EEGNet accuracy for all within session (black), cross session (red), cross
subject (blue) cases across several amounts of sensor failure (0%, 10%, 30%,
70%). Ovals in 0% failure case represent standard error of mean (SEM) across 5
model random seeds. Ovals in 10%, 30%, 70% sensor failure plots represent SEM
across 3 sensor failure random seeds. (b) Decoding accuracy of TOTEM (solid
line) and EEGNet (dashed line) when trained on B1, and tested against within
session (black), cross session (red), and cross subject (blue) cases across
0-100% sensor failure percentages. Additional performances can be found in
Appendix 5.1.

When analyzing the 0% sensor failure generalizability cases, we see TOTEM
outperforms EEGNet on cross session and subject cases, and match EEGNet on
within session performance. On a case by case basis, we see that TOTEM beats
EEGNet on all cross session modeling experiments, and beats in 4/8 cross
subject modeling experiments. The within session experiments all cluster
around the equal performance line, which means that TOTEM is able to do just
as well as EEGNet when evaluated at the baseline case. The difference in
generalizability performance between the two models while within session
decoding is fairly consistent shows that it is important to evaluate models
along generalization dimensions when selecting for models.

When analyzing the 10%, 30%, and 70% generalizability cases, we see that for
nearly all modeling experiments, the same generalizbility results we saw hold,
with TOTEM’s performance becoming relatively better than EEGNet’s as more
sensors are failed. This can be observed by the dots shifting leftward and
more into the green regions as more sensors are failed (Figure 5). The within
session decoding accuracies also favor TOTEM as more sensors are failed, which
is beneficial for modeling within session datasets where sensor failure may
appear randomly within a session.

When analyzing the trends of each model against a finer resolution of sensor
failure percentages, we see that TOTEM maintains a higher decoding accuracy
for longer, while EEGNet has a more linear declining performance with
additional sensors removed (Figure 5 (b)). This is observed even when both
models start with the same decoding accuracy, so it is not due to TOTEM simply
generalizing better in cross session and cross subject cases. Additional
results trained on each of our 4 sessions can be found in Appendix 5.1.

These results suggest that (1) TOTEM’s tokenization + transformers approach
create more generalizable representations when compared to EEGNet’s CNN
kernels, and (2) studying generalization highlights model performance
differences as the within session conclusions differ from the generalizability
conclusions.

###  3.2 Latent Codebook Analysis

Figure 6: Codebook similarity matrices and visualization of codewords within
codebooks. (a) Matrices of MSE between matching codewords in the two labeled
codebooks. Codewords are ordered roughly by within session mean MSE. Colored
arrows along axis denote visualized codewords highlighted in (b). (b) All
codewords in each codebook. Highlighted are three examples of matched
codewords between codebooks. MSE values between codewords are: purple=0.065,
pink=0.018, orange=0.013.

To better understand TOTEM’s generalization capabilities, we investigate its
learned latent codebook across sessions, subjects, and sensor availability. If
the codebooks learned are generalizable, we would expect that the codewords
have corresponding matches to the codewords in the other codebooks, allowing a
new recording session to be reconstructed and represented with high fidelity.
Indeed, when we plot the MSE between codewords across different codebooks
(Figure 6(a) Cross Subject), we find that the codebook similarity matrices are
qualitatively comparable to within codebook (Figure 6(a), Within Session),
especially when we match each codeword from the new codebook (A2) with its
minimum MSE codeword in the original codebook (B1). The matched codewords also
look highly similar across codebooks (Figure 6(b)) and are quantitatively low
in MSE. Across the whole codebook, the average MSE of matched codewords is
near zero and a majority of codewords are used when matching (Table 1),
suggesting high overlap between codebooks learned across generalizability
cases. All generalization cases studied here learn similar codebooks according
to these qualitative and quantitative metrics (additional visualizations in
Appendix 5.2). This shows that TOTEM’s learned tokenization space may be able
to learn a generalized codebook for our time series data modality of interest,
potentially allowing for mass training of the downstream transformer encoders
using the same tokenized latent space.

| Average MSE | Num. Subselected Codewords
---|---|---
Baseline Case |  |
Within Session: B1 $\rightarrow$ B1 | 0.0 | 256
Generalization Cases |  |
Cross Session: B$1_{sub}$ $\rightarrow$ B2 | 0.06 | 163
Cross Subject: B$1_{sub}$ $\rightarrow$ A2 | 0.05 | 188
10% sensors: B$1_{sub}$ $\rightarrow$ B1 10% | 0.04 | 194
Table 1: Average MSE of matched codewords in baseline case and generalization
cases. Number of Subselected Codewords denote the number of codewords from the
original codebook (B1) that were matched to the generalization codebook.

##  4 Conclusion + Limitations + Future Work

We find that tokenization + transformers are a promising approach to modeling
time series neural data which have high variability in datasets and
recordings. Specifically, we demonstrate that compared with one of the most
performant and popular CNN models, tokenization+transformers outperform in
numerous generalization cases. These models are also ripe for interpretability
analysis which can uncover new findings about time series neural data. Our
study further shows how important it is to consider generalizability cases
when selecting models.

Currently we only model sensor failure as zeroing out data. This masking
technique is widely used machine learning, but does not cover all sensor
failure cases. Creating more sensor failure modes is a meaningful direction
for future work. Extending our interpretability analysis to the downstream
transformer is another valuable future direction. We would also like to adapt
the framework to work with more types of neural time series data such as
sparsely sampled neural time series data such as LFP from neuropixels, Utah
Arrays, and stereoencepholography (sEEG) recordings, and more varied layout
modalities such as electrocortiography (ECoG). This work could enable
automatic noisy sensor detection, interesting interpretations of temporal and
spatial dimensions, and building of foundation models for neural time series.

## References

  * Aaron et al. (2017) Aaron, Oriol Vinyals, and Koray Kavukcuoglu.  Neural discrete representation learning, 2017.  URL https://arxiv.org/abs/1711.00937v2.
  * Ang et al. (2012) Kai Keng Ang, Zheng Yang Chin, Chuanchu Wang, Cuntai Guan, and Haihong Zhang.  Filter bank common spatial pattern algorithm on bci competition iv datasets 2a and 2b.  _Frontiers in Neuroscience_ , 6, Jan 2012.  doi: https://doi.org/10.3389/fnins.2012.00039.  URL https://www.frontiersin.org/journals/neuroscience/articles/10.3389/fnins.2012.00039/full.
  * García-Murillo et al. (2023) Daniel Guillermo García-Murillo, Andrés Marino Álvarez Meza, and Cesar German Castellanos-Dominguez.  Kcs-fcnet: Kernel cross-spectral functional connectivity network for eeg-based motor imagery classification.  _Diagnostics_ , 13(6):1122–1122, Mar 2023.  doi: https://doi.org/10.3390/diagnostics13061122.  URL https://www.mdpi.com/2075-4418/13/6/1122.
  * Gonschorek et al. (2021) Dominic Gonschorek, Larissa Höfling, Klaudia P Szatko, Katrin Franke, Timm Schubert, Benjamin Dunn, Philipp Berens, David Klindt, and Thomas Euler.  Removing inter-experimental variability from functional data in systems neuroscience.  _Advances in Neural Information Processing Systems_ , 34:3706–3719, Dec 2021.  URL https://proceedings.neurips.cc/paper_files/paper/2021/hash/1e5eeb40a3fce716b244599862fd2200-Abstract.html.
  * Gwon et al. (2023) Daeun Gwon, Kyungho Won, Minseok Song, Chang S Nam, Sung Chan Jun, and Minkyu Ahn.  Review of public motor imagery and execution datasets in brain-computer interfaces.  _Frontiers in Human Neuroscience_ , 17, Mar 2023.  doi: https://doi.org/10.3389/fnhum.2023.1134869.  URL https://www.frontiersin.org/articles/10.3389/fnhum.2023.1134869/full#B76.
  * Kaya et al. (2018) Murat Kaya, Mustafa Kemal Binli, Erkan Ozbay, Hilmi Yanar, and Yuriy Mishchenko.  A large electroencephalographic motor imagery dataset for electroencephalographic brain computer interfaces.  _Scientific Data_ , 5(1), Oct 2018.  doi: https://doi.org/10.1038/sdata.2018.211.  URL https://www.nature.com/articles/sdata2018211.
  * Krumpe et al. (2017) T Krumpe, K Baumgärtner, W Rosenstiel, and M Spüler.  Non-stationarity and inter-subject variability of eeg characteristics in the context of bci development.  _7th Graz Brain-Computer Interface Conference 2017, Graz, Austria_ , 2017.  doi: http://hdl.handle.net/10900/83769.  URL https://tobias-lib.ub.uni-tuebingen.de/xmlui/handle/10900/83769.
  * Lawhern et al. (2018) Vernon J Lawhern, Amelia J Solon, Nicholas R Waytowich, Stephen M Gordon, Chou P Hung, and Brent J Lance.  Eegnet: a compact convolutional neural network for eeg-based brain–computer interfaces.  _Journal of Neural Engineering_ , 15(5):056013, jul 2018.  doi: 10.1088/1741-2552/aace8c.  URL https://dx.doi.org/10.1088/1741-2552/aace8c.
  * Parvizi & Kastner (2018) Josef Parvizi and Sabine Kastner.  Promises and limitations of human intracranial electroencephalography.  _Nature Neuroscience_ , 21(4):474–483, Mar 2018.  doi: https://doi.org/10.1038/s41593-018-0108-2.  URL https://www.nature.com/articles/s41593-018-0108-2.
  * Peterson et al. (2021) Steven M Peterson, Zoe Steine-Hanson, Nathan Davis, Rajesh P N Rao, and Bingni W Brunton.  Generalized neural decoders for transfer learning across participants and recording modalities.  _Journal of Neural Engineering_ , 18(2):026014, mar 2021.  doi: 10.1088/1741-2552/abda0b.  URL https://dx.doi.org/10.1088/1741-2552/abda0b.
  * Saha & Baumert (2020) Simanto Saha and Mathias Baumert.  Intra- and inter-subject variability in eeg-based sensorimotor brain computer interface: A review.  _Frontiers in Computational Neuroscience_ , 13, Jan 2020.  doi: https://doi.org/10.3389/fncom.2019.00087.  URL https://www.frontiersin.org/articles/10.3389/fncom.2019.00087/full.
  * Stieger et al. (2021) James R Stieger, Stephen A Engel, and Bin He.  Continuous sensorimotor rhythm based brain computer interface learning in a large population.  _Scientific Data_ , 8(1), Apr 2021.  doi: https://doi.org/10.1038/s41597-021-00883-1.  URL https://www.nature.com/articles/s41597-021-00883-1.
  * Talukder et al. (2022) Sabera Talukder, Jennifer J. Sun, Matthew Leonard, Bingni W. Brunton, and Yisong Yue.  Deep neural imputation: A framework for recovering incomplete brain recordings, 2022.  URL https://doi.org/10.48550/arXiv.2206.08094.
  * Talukder et al. (2024) Sabera Talukder, Yisong Yue, and Georgia Gkioxari.  Totem: Tokenized time series embeddings for general time series analysis, 2024.  URL https://arxiv.org/abs/2402.16412v1.
  * Tangermann et al. (2012) Michael Tangermann, Klaus Robert Müller, Ad Aertsen, Niels Birbaumer, Christoph Braun, Clemens Brunner, Robert Leeb, Carsten Mehring, Kai J Miller, Gernot R Müller-Putz, Guido Nolte, Gert Pfurtscheller, Hubert Preissl, Gerwin Schalk, Alois Schlögl, Carmen Vidaurre, Stephan Waldert, and Benjamin Blankertz.  Review of the bci competition iv.  _Frontiers in Neuroscience_ , 6, Jan 2012.  doi: https://doi.org/10.3389/fnins.2012.00055.  URL https://www.frontiersin.org/journals/neuroscience/articles/10.3389/fnins.2012.00055/full?ref=https%3A%2F%2Fgithubhelp.com.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin.  Attention is all you need, 2017.  URL https://arxiv.org/abs/1706.03762.
  * Xu et al. (2020) Lichao Xu, Minpeng Xu, Yufeng Ke, Xingwei An, Shuang Liu, and Dong Ming.  Cross-dataset variability problem in eeg decoding with deep learning.  _Frontiers in Human Neuroscience_ , 14, Apr 2020.  doi: https://doi.org/10.3389/fnhum.2020.00103.  URL https://www.frontiersin.org/articles/10.3389/fnhum.2020.00103/full.
  * Xu et al. (2021) Lichao Xu, Minpeng Xu, Zheng Ma, Kun Wang, Tzyy-Ping Jung, and Dong Ming.  Enhancing transfer performance across datasets for brain-computer interfaces using a combination of alignment strategies and adaptive batch normalization.  _Journal of Neural Engineering_ , 18(4):0460e5–0460e5, Aug 2021.  doi: https://doi.org/10.1088/1741-2552/ac1ed2.  URL https://iopscience.iop.org/article/10.1088/1741-2552/ac1ed2.

##  5 Appendix

###  5.1 Sensor failure performance across subjects

Figure 7: Classifier performance under sensor failure. A) Classifier accuracy
across variable sensor failure %s. Each column represents models trained on
the dataset in the title. Error bars represent SEM of 5 random seeds which
randomly select the failed sensors. wSes=Within Session, xSes=Cross Session,
xSub=Cross Subject. B) TOTEM & EEGNet accuracy difference with mean of
differences in red. If a point lies above 0, TOTEM decodes better.

###  5.2 Latent codebook generalizability across all cases

(a) Within Session Codebook MSE Matrices

(b) Matrices of MSE between matching codewords in the two labeled codebooks.
Codewords are ordered roughly by within session mean MSE. Colored arrows along
axis denote highlighted codewords in (c).

(c) All codewords in each codebook. Highlighted are three examples of matched
codewords between codebooks

Figure 8: Additional latent codebook visualizations


# Finite element method modeling of expansion of irradiated rocks: focusing on
the minerals

Y. Oida  MRI Research Associates, Inc., Chiyoda, Tokyo 100-0014 Japan.  S.
Sakuragi<EMAIL_ADDRESS>MRI Research Associates, Inc.,
Chiyoda, Tokyo 100-0014 Japan.  T. Igari  MRI Research Associates, Inc.,
Chiyoda, Tokyo 100-0014 Japan.  Y. Nakajima  MRI Research Associates, Inc.,
Chiyoda, Tokyo 100-0014 Japan.  E. Ro  MRI Research Associates, Inc., Chiyoda,
Tokyo 100-0014 Japan.  F. Shimizu  MRI Research Associates, Inc., Chiyoda,
Tokyo 100-0014 Japan.  M. Futami  MRI Research Associates, Inc., Chiyoda,
Tokyo 100-0014 Japan.  Y. Hakozaki  Graduate School of Engineering, Chiba
University, Inage, Chiba 263-8522, Japan.  T. Ohkubo  Graduate School of
Engineering, Chiba University, Inage, Chiba 263-8522, Japan.  H. Ishikawa
Department of Mechanical Systems Engineering, Tokyo University of Agriculture
and Technology, Koganei, Tokyo 184–8588, Japan.  S. Takada  Department of
Mechanical Systems Engineering, Tokyo University of Agriculture and
Technology, Koganei, Tokyo 184–8588, Japan.  Institute of Engineering, Tokyo
University of Agriculture and Technology, Koganei, Tokyo 184–8588, Japan.  S.
Sawada  Kajima Corporation, Minato, Tokyo 107-8388 Japan.  K. Murakami
Graduate School of Engineering, The University of Tokyo, Bunkyo, Tokyo
113-8656, Japan.  K. Suzuki  Mitsubishi Research Institute, Inc., Chiyoda,
Tokyo 100-0014 Japan.  I. Maruyama<EMAIL_ADDRESS>Graduate School of Engineering, The University of Tokyo, Bunkyo, Tokyo
113-8656, Japan.  Graduate School of Environmental Studies, Nagoya University,
Chikusa, Nagoya 464-8603, Japan.

###### Abstract

Finite element method (FEM) modeling of the volumetric expansion phenomenon
associated with the accumulation of irradiation was performed on rocks in a
concrete for nuclear power plant. The FEM mesh of sandstone, tuff, and granite
was created based on microscopic images, and the volumetric expansion due to
irradiation was calculated. The saturated value of the volume expansion due to
irradiation accumulation was calculated, and it was shown that the saturated
value of the volume expansion was well explained by the experimental value. In
addition, the FEM analysis indicates that irradiation-induced degradation of
rock propagates through a localized concentration of stress followed by
fracture with cracks spreading throughout the rock.

##  I Introduction

In Japan, nuclear power plants have been in continuous operation for 30 years,
and there is a need for a method to evaluate the deterioration of power plant
buildings from the viewpoint of safety. Since nuclear power plant buildings
are concrete structures, the safety of the buildings will be determined by
evaluating the deterioration of the concrete, as is the case with ordinary
buildings. On the other hand, concrete used in nuclear power plants is under a
special environment where it is constantly exposed to neutrons for several
decades, and thus, the degradation phenomena of general structures cannot be
appropriated and discussed in the same way. Therefore, there is a need to
clarify the mechanism of degradation phenomena of concrete caused by neutron
irradiation in order to obtain general knowledge on the safety of nuclear
power plants [1, 2, 3].

In previous reports, irradiation degradation of concrete has been discussed
from an experimental point of view: expansion and degradation due to
irradiation accumulation in the rocks configuring the concrete [4, 5, 6].
Dimensional measurements have revealed the existence of irradiation-induced
expansion from a macroscopic viewpoint, and it has been reported that volume
expansion of up to 20% or more occur in sandstone, tuff, and granite [7]. In
addition, x-ray diffraction (XRD) measurement on rocks have also discussed
degradation phenomena from a microscopic viewpoint, such as changes in the
lattice constants of mineral crystals constituting the rocks, and they reveals
that these irradiation expansions are caused by the amorphization of mineral
crystals associated with neutron irradiation [7, 8, 9, 10]. On the other hand,
even for quartz with high radiosensitivity, the volume expansion obtained from
XRD is about $\sim$18%, which cannot explain the macroscopic results based on
dimensional measurements. This phenomenon is qualitatively understood to be
the expansion of the mineral crystals constituting the rocks, which causes
cracks throughout the rocks and thus expands the rocks beyond the effect of
amorphization. However verifying this phenomenon experimentally, experiments
using neutron reactors have not provided systematic experimental data due to
the difficulty of conducting systematic experiments on a large number of rock
samples, and thereby a universal mechanism explaining the irradiation swelling
of rock from the mineral perspective has not been clarified.

In contrast to clarifying the universal curve of irradiation degradation in
rocks depending on the constituent mineral species, an approach utilizing
computer science is useful for systematically discussing the relationship
between the expansion of various minerals and rock degradation. Therefore, as
a preparatory step, this study examined a method to create an FEM model
focusing on the minerals constituting the rock from photographs of the actual
rock. The FEM model considering the fracture of grain boundaries for
sandstone, tuff, and granite reproduced well the maximum value of irradiation
expansion obtained from the experiment. We propose that the FEM model can be
used to understand the irradiation degradation and expansion of rockss using
computational methods.

##  II modelling and method

Ansys LS-DYNA, a nonlinear structural analysis software developed by Livermore
Software Technology Corporation, was used for the FEM analysis [11]. Since
rock fracture is well explained by Griffith theory [12, 13], the FEM
constitutive equation (spring-damper model: MAT189 in LS-DYNA) was used, and
modeling was performed considering fracture at the grain boundary. The
expansion process associated with the accumulation of irradiation was treated
as thermal expansion reaching thermal equilibrium and calculated by FEM. FEM
mesh elements that were determined to be destructive were deleted to express
the destruction.

Polarized light microscope images were taken of sandstone, tuff, and granite
surfaces to create the FEM mesh. The mineral species were identified by image
analysis using color resolved [14] watershed algorithm, and voxel meshes were
generated {Fig. 1 (a) and (b)}. The percentages of each mineral in the rock
obtained from image analysis are shown in Table 1. Note that this value agrees
well with the value obtained from the XRD measurement within an error of
$\sim$1%. For the irradiation dependence of coefficient of linear expansion,
values reported in Denisov’s report were used [7]. The physical properties of
the minerals were obtained from molecular dynamics (MD) calculations [15, 16].
Since the crystallographic orientation of the minerals cannot be determined
from polarized light microscope images, the coefficients of linear expansion
were set by randomly assigning crystallographic orientations to the minerals
in this model.

And then, fracture stresses were determined so that the FEM Brazillian test
reproduced the tensile strength values obtained from reports of previous
experiments. For FEM analysis of the Brazilian test, a circularly cut mesh of
polarized light microscope images [17] was used. The lower end of the circle
was constrained and forced displacement was applied from the upper end to
simulate the Brazilian test. For simplicity, the vertical and shear failure
stresses were assumed to be the same in this model. Although the maximum
friction coefficient in rocks is generally known to be around 0.6 [18, 19], an
overestimation of friction was avoided by uniformly using a friction
coefficient = 0.3. Rock expansion due to irradiation degradation was evaluated
by performing irradiation expansion calculations using the above parameters
for model.

Figure 1:  (a) Polarization microscope image of sandstone. (b) Results of mineral and grain boundary identification by image analysis. Pink indicates quartz, yellow indicates feldspar, and dark blue indicates biotite.  Table 1:  Percentage of minerals composing the rock.  | Sandstone | Tuff | Granite
---|---|---|---
Quartz | 58.3 % | 88.2 % | 34.1 %
Feldspar | 23.4 % | - | 57.6 %
Chlorite | 18.3 % | 3.9 % | -
Biotite | - | 7.9 % | 8.3 %

##  III Result and Discussion

Table 2 shows the fracture stress parameters in the FEM that well explain the
experimental values of tensile strength in the sandstone and tuff, with errors
of 6% and 1%, respectively [20, 21]. The tensile strength of the granite was
not reproduced by the analysis even by using the reported experimental values
of fracture stress, and the tensile strength in the FEM model showed excessive
values even when the friction coefficient was set to as small as 0.001.

Table 2:  FEM parameters of rocks obtained from Brazilian tests.  | Sandstone | Tuff | Granite
---|---|---|---
Fracture stress (MPa) | 12 | 20 | 5
Friction coefficient | 0.3 | 0.3 | 0.001
Reaction force (kN) | 0.0031 | 0.0033 | 0.00915
Tensile strength (MPa) | 19.1 | 20.3 | 56.3

To investigate the tendency of failure in the Brazilian test, the principal
stresses for each rock were examined when the reaction force reached its
maximum value in the Brazilian test. Principal stress contour plots for each
rock are shown in Fig. 2. In the tuff with the smallest error in tensile
strength, the stresses are linearly generated from the point of forced
displacement to the point of constraint. On the other hand, in granite, the
stress is anisotropic and excessive stress is generated in the mineral grains
of the granite. This excess stress is expected to be resolved by fracture
within the mineral grains of the granite, which is not originally assumed in
this FEM model. Alternatively, it could be resolved by performing the analysis
for even larger granites and creating a model in which a grain boundary
appears between the forced displacement point and the constraint point.
However, the actual size of the mineral grains in the experiment is unknown,
and it is also unknown whether fracture occurred inside the mineral grains in
the experiment. The results of this FEM analysis indicate that the
relationship between rock grain size and tensile strength should be carefully
investigated in experiments.

Figure 2:  Contour plot of maximum stress during Brazilian test.

We treated irradiation accumulation as a thermal expansion phenomenon in a
sufficiently thermal equilibrium state, and analyzed the irradiation expansion
phenomenon by FEM using the above parameters. The volumetric expansion rate
$\Delta V$ and the linear strain $\varepsilon$ of the rock resulting from
irradiation accumulation are defined as follows:

$\displaystyle\Delta V$
$\displaystyle=\frac{3}{2}\\{\left(1+\varepsilon\right)^{2}-1\\}+1,\varepsilon$
$\displaystyle=\frac{{\mathrm{\Delta
l}}_{ave}}{\sqrt{\left(\frac{X}{2}\right)^{2}+\left(\frac{Y}{2}\right)^{2}}},$
(1)

where ${\mathrm{\Delta l}}_{ave}$ is the average displacement of the nodes at
each vertex of the FEM mesh, X is the length of the analytical model in the
$x-$direction, and Y is the length in the $y-$direction. The maximum, minimum,
and average values of the volumetric expansion of each rock for each neutron
fluence are calculated and are shown in Fig 3. Experimental values of the
irradiation expansion for sandstone and granite have been reported by Denisov
et al. [7], and our FEM calculations reproduce the absolute values and
qualitative trends well.

According to Denisov et al.[7], for quartz-rich rocks in general, the trend of
irradiation expansion can be expressed by the following:

$\displaystyle\frac{\mathrm{\Delta V}}{V}=\frac{a{\frac{\mathrm{\Delta
V}}{V}}_{max}\left(e^{bn_{cm}}-1\right)}{{\frac{\mathrm{\Delta
V}}{V}}_{max}+ae^{bn_{cm}}},$ (2)

where $a$ and $b$ are parameters that depend on the mineral species,
$\left(\mathrm{\Delta V}/V\right)_{max}$ is the maximum expansion rate of the
mineral, and $n_{cm}$ is the number of displaced atoms, proportional to the
irradiation dose. This function shows an S-shaped curve, which is consistent
with the irradiation expansion results shown in the Fig. 3.

Figure 3:  Results of FEM analysis of volumetric expansion rate with
irradiation accumulation in rocks. The solid line indicates the mean value,
and the dashed lines indicate the maximum and minimum values, respectively.

The maximum volumetric expansion of individual minerals calculated from the
linear expansion coefficients used in our input is 17.9% for quartz [7]. The
maximum expansion is higher than that of quartz in all the rocks analyzed in
this study. These results suggest that the volumetric expansion of rocks due
to radiation accumulation is not only caused by linear expansion of minerals,
but also by fracture at grain boundaries between minerals, which creates voids
that increase the volumetric expansion of rocks. In fact, the present analysis
shows that fracture occurs in the FEM mesh at an irradiation dose of $1\times
10^{24}neutron/m^{2}$. In addition, the maximum expansion of the rock tends to
be higher for rocks with higher quartz content in the range analyzed in this
study.

In order to understand the phenomenon of irradiation expansion in rocks,
stress contour plots were output for each neutron fluence. The stress contour
diagram for tuff is shown in Fig. 4. Note that, stress contour images for
sandstone and granite are similar to those for tuff. From the stress contour
diagram, it can be seen that at neutron fluence $0.5\times
10^{24}neutron/m^{2}$, a localized spot of fracture stress is observed in the
interior of the rock, which tends to spread over the entire rock from fluence
$1-3\times 10^{24}neutron/m^{2}$. The irradiation doses from $1-3\times
10^{24}neutron/m^{2}$ correspond to an increase in the linear expansion
coefficient of quartz. Although the coefficient of linear expansion of quartz
decreases and saturates after $3\times 10^{24}neutron/m^{2}$, the expansion of
the rock continues to increase. This may be attributed to the fact that the
fracture of the rock allows voids to be generated at the grain boundary, and
these voids gradually develop. Based on the above results, the following
three-step model can be considered as the irradiation expansion mechanism of
the rock: (1) $\sim 0.5\times 10^{24}neutron/m^{2}$: irradiation expansion due
to the superposition of expansion of each mineral that does not cause fracture
inside the rock, (2) $1-3\times 10^{24}neutron/m^{2}$: irradiation expansion
with fracture inside the rock due to rapid expansion of quartz, (3) over
$3\times 10^{24}neutron/m^{2}$: irradiation swelling due to the generation of
cracks. The result of our FEM analysis gives a reasonable interpretation of
the curve drawn by functional of Eq. 2.

Figure 4:  Dependence of irradiation accumulation on stresses inside tuff
obtained from FEM analysis.

Our present FEM model does not take into account the fracture inside the
mineral, as mentioned above. The fracture of the mineral interior is expected
to affect the generation and relaxation of local stress, especially in low-
irradiation regions, as shown in the stress contour diagram (Fig. 4).
Therefore, with further detailed experiments, a more accurate FEM model will
be considered by modeling the fracture phenomena associated with irradiation
accumulation based on the experimental results, especially in the low-
irradiation regions.

##  IV conclusion

Our present FEM modeling shows that it is possible to reproduce the fracture
and expansion phenomena associated with irradiation accumulation of rocks used
in concrete buildings in large radiation accumulation regions through a two-
dimensional FEM analysis using experimental values of linear expansion
coefficients of minerals as inputs. Using this model, it will be possible to
clarify radiation-induced degradation phenomena for a wide variety of rocks
from the viewpoint of computer science. The present analysis also indicates
that careful investigation of fracture phenomena inside minerals is necessary,
especially in mechanical tests and irradiation experiments in low-irradiation
regions. Clarifying the fracture phenomena in regions of small irradiation
accumulation in rocks are important from the viewpoint of maintenance of
nuclear power plant. It is expected that our model will be refined by feedback
from experiments, and that the safety of nuclear power plants will be
discussed from computer simulations in the future.

## acknowledgments

This work constitutes part of the Japan Concrete Aging Management Program on
Irradiation Effects (JCAMP), sponsored by METI in Japan.

## References

  * Braverman  _et al._ [2004] J. Braverman, C. Miller, C. Hofmayer, B. Ellingwood, D. Naus, and T. Chang, Degradation assessment of structures and passive components at nuclear power plants, Nuclear engineering and design 228, 283 (2004).
  * Le Pape  _et al._ [2020] Y. Le Pape, J. Sanahuja, and M. H. Alsaid, Irradiation-induced damage in concrete-forming aggregates: revisiting literature data through micromechanics, Materials and Structures 53, 1 (2020).
  * Maruyama  _et al._ [2021] I. Maruyama, M. Takizawa, O. Sato, J. Etoh, O. Kontani, S. Sawada, and S. Ishikawa, Post-japanese nra research project on soundness evaluation criteria for radiation-induced concrete degradation, in  _Nuclear Power Plant Life Management. Proceedings of an International Conference. Supplementary Files_ (2021).
  * Hilsdorf  _et al._ [1978] H. Hilsdorf, J. Kropp, and H. Koch, The effects of nuclear radiation on the mechanical properties of concrete, Special Publication 55, 223 (1978).
  * Maruyama  _et al._ [2018] I. Maruyama, S. Ishikawa, J. Yasukouchi, S. Sawada, R. Kurihara, M. Takizawa, and O. Kontani, Impact of gamma-ray irradiation on hardened white portland cement pastes exposed to atmosphere, Cement and Concrete Research 108, 59 (2018).
  * Pape  _et al._ [2018] Y. L. Pape, M. H. F. Alsaid, and A. B. Giorla, Rock-forming minerals radiation-induced volumetric expansion – revisiting literature data, Journal of Advanced Concrete Technology 16, 191 (2018).
  * Denisov  _et al._ [2012] A. V. Denisov, V. B. Dubrovskii, and V. N. Solovyov, Radiation resistance of mineral and polymer construction materials, ZAO MEI Publishing House (in Russian) (2012).
  * Devine [1994] R. Devine, Macroscopic and microscopic effects of radiation in amorphous sio2, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 91, 378 (1994).
  * Devine [1990] R. Devine, On the physical models of annealing of radiation induced defects in amorphous sio2, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 46, 261 (1990).
  * Hobbs [1995] L. W. Hobbs, The role of topology and geometry in the irradiation-induced amorphization of network structures, Journal of non-crystalline solids 182, 27 (1995).
  * [11] LS-DYNA, https://www.lstc.com/products/ls-dyna.
  * Griffith [1921] A. A. Griffith, Vi. the phenomena of rupture and flow in solids, Philosophical transactions of the royal society of london. Series A, containing papers of a mathematical or physical character 221, 163 (1921).
  * Liebowitz and Sih [1968] H. Liebowitz and G. Sih,  _Mathematical theories of brittle fracture._ (ACADEMIC PRESS, INC., 1968).
  * Chen  _et al._ [2004] S. Chen, Z. Yue, and L. Tham, Digital image-based numerical modeling method for prediction of inhomogeneous rock failure, International journal of rock mechanics and mining sciences 41, 939 (2004).
  * Okada  _et al._ [2020] N. Okada, T. Ohkubo, I. Maruyama, K. Murakami, and K. Suzuki, Characterization of irradiation-induced novel voids in $\alpha$-quartz, AIP Advances 10, 125212 (2020).
  * [16] to be submitted.
  * [17] private communication.
  * Tullis [2007] T. Tullis, Friction of rock at earthquake slip rates, Earthquake Seismology 4, 131 (2007).
  * Byerlee [1978] J. Byerlee, Friction of rocks, pure and applied geophysics 116, 615 (1978).
  * N.Takashiki [1987] N.Takashiki, On the deformation and failure behaviour of omotani rhyolite (in japanese), Japan Society of Civil Engineers 19, 31 (1987).
  * Izawa  _et al._ [2020] S. Izawa, S. Ogata, H. Yasuhara, N. Kinoshita, and K. Kishida, Evaluation of fracture evolution of granite during brazilian test by numerical analysis of fracturing process considering mineral distribution, Journal of the Society of Materials Science Japan 69, 236 (2020).