Modeling heterogeneity in higher-order moments
while preserving mean and variance:
application to spatio-temporal modeling

The modeling of higher-order moments, such as skewness and kurtosis, is crucial in many real-world applications (Tagle et al., 2019). For example, infection risk may exhibit a skewed distribution with extremely high values during epidemic periods, and the strength of skewness can vary across different regions and time periods. Therefore, describing the process using a non-Gaussian model that can flexibly handle such heterogeneous skewness and higher-order moments is desirable.

To account for non-Gaussianity, a wide variety of methods have been developed, as reviewed by Yan et al. (2020). A prominent example is the closed-skew normal (CSN) distribution, known for its expressive power (Azzalini, 1985; Azzalini and Capitanio, 1999). However, replacing the normal distribution in a model with a CSN distribution often disrupts the mean and variance properties, which are essential for interpretability. In particular, spatio-temporal models—actively studied in fields such as disease mapping (Best et al., 2005), species distribution modeling (Dormann et al., 2007), and socio-economic analysis (Anselin, 2022) —rely on the mean regression function to interpret the influence of the specified covariates and covariance specified to quantify the strength of spatial and temporal correlations.

Márquez-Urbina and González-Farías (2022) addressed this issue by proposing a flexible subclass of the CSN distribution that maintains the mean and variance. This property enhances the interpretability of parameters such that when the mean is specified by a regression function, the coefficients and variance parameters can be interpreted similarly to those in a conventional Gaussian regression model. Additionally, Tan and Chen (2024) proposed an extended model known as the closed skew normal subclass (CSNS) distribution, which addresses the heterogeneity of higher-order moments. The CSNS distribution still has limitations in the range of achievable skewness and kurtosis, restricting its representational flexibility. Moreover, Tan and Chen (2024) used the CSNS distribution as an approximation in variational inference, rather than as a component of modeling. Notably, although the model’s properties vary depending on the matrix decomposition method used internally, it is restricted to LU decomposition. Similarly, no established Bayesian estimation method is available despite its potential merits for modeling uncertainty in spatio-temporal phenomena (Cressie and Wikle, 2011; Márquez-Urbina and González-Farías, 2022).

To overcome these limitations, we propose a general model that encompasses the CSNS distribution as a special case and accommodates the heterogeneity of higher-order moments while preserving the mean and variance. Subsequently, several properties of the proposed model are presented. In particular, since the characteristics of the model vary depending on the matrix factorization method employed, we explain the specific matrix factorization techniques suitable for adapting to spatio-temporal data. Furthermore, we develop an efficient Bayesian inference method that is particularly effective for temporal data.

The remainder of this paper is organized as follows. Section 2 provides background on the CSN and CSNS distributions. Section 3 introduces the proposed model and its theoretical properties. Section 4 investigates the impact of different matrix factorization methods, particularly in spatio-temporal context. Section 5 develops an efficient Bayesian inference method. Section 6 presents the simulations conducted to the model estimation accuracy, and Section 7 applies the model to the estimation of production functions across U.S. states.

The CSN distribution is an extension of the multivariate normal distribution that incorporates higher-order moments. For $\bm{s}\sim\mathcal{CSN}_{N,M}(\bm{\mu},\bm{\Sigma},\bm{D},\bm{\nu},\bm{\Omega})$, the probability density function is given by:

Here, $\bm{\mu}\in\mathbb{R}^{p}$, $\bm{\nu}\in\mathbb{R}^{q},\bm{D}\in\mathbb{R}^{q\times p}$, and $\bm{\Sigma}\in\mathbb{R}^{N\times N},\bm{\Omega}\in\mathbb{R}^{M\times M}$ are covariance matrices. $\phi_{N}(\bm{s};\bm{\mu},\bm{\Sigma})$ and $\Phi_{N}(\bm{s};\bm{\mu},\bm{\Sigma})$ denotes the probability distribution function and cumulative distribution function of an $N$-dimensional normal distribution with mean $\bm{\mu}$ and covariance $\bm{\Sigma}$, respectively. Finally, $\bm{0}_{M}$ represents an $M$-dimensional zero vector.

As a modified version of the CSN distribution, the CSNS distribution ensures that the mean and variance remain consistent with those of a multivariate normal distribution. The CSNS distribution is derived by first defining the following random variable:

where $\mathcal{N}_{N}(\bm{\mu},\bm{\Sigma})$ denotes a multivariate normal distribution with mean $\bm{\mu}$ and covariance $\bm{\Sigma}$, and $\mathcal{TN}(\mu,\sigma^{2})$ denotes a normal distribution $\mathcal{N}_{1}(\mu,\sigma^{2})$ truncated below $0$. The mean of the above random variable is $\sqrt{2/\pi}\theta_{i}$ and the variance is $f_{i}^{-2}=(\pi+(\pi-2)\theta_{i}^{2})/\pi$. To standardize the variable, we define $t_{i}=f_{i}(\theta_{i}(u_{i}-\sqrt{2/\pi})+v_{i})$, ensuring that $t_{i}$ has a mean of $0$ and a variance of $1$. Setting

guarantees that $\bm{s}$ maintains mean $\bm{\mu}$ and covariance $\bm{\Sigma}$. Here, $\bm{s}=(s_{1},\ldots,s_{N})^{\top}$, $\bm{t}=(t_{1},\ldots,t_{N})^{\top}$, and $\bm{\Sigma}^{LU}$ is the LU decomposition of $\bm{\Sigma}$. The distribution of $\bm{s}$ can be expressed as the CSN distribution:

where $\bm{\mu}_{\bm{u}}=(\sqrt{2/\pi}\theta_{1},\ldots,\sqrt{2/\pi}\theta_{N})^{\top}$, $\bm{F}=diag(f_{1},\ldots,f_{N})$, $\bm{\Theta}=diag(\theta_{1},\ldots,\theta_{N})$, and $\bm{I}_{N}$ is the identity matrix.

Figure 1 illustrates the probability density function of the CSNS distribution with mean $0$ and variance $1$ for different values of $\theta$. When $\theta=0$, the CSNS distribution coincides with the normal distribution. As $\theta$ increases, the skewness increases, and the distribution approaches a shifted half-normal distribution.

In this section, a general model that can handle heterogeneity in higher-order moments while preserving the mean and variance is proposed, based on the idea of the CSNS distribution. In Section 3.1, we provide an overview of the model along with several illustrative examples. Subsequently, in Section 3.2, we discuss the model’s closure under linear transformations as well as skewness and kurtosis.

In the proposed model, characteristics vary based on the matrix decomposition method. Although the scope of this model is not necessarily limited to spatial or temporal data, it is well-suited for explaining changes in model properties. Additionally, since spatial and temporal models are often designed based on the mean and variance, preserving these statistical properties. In Section 4.1, we consider applications to spatial models, followed by applications to temporal models in Section 4.2.

This section proposes a Bayesian estimation method for the parameters of the proposed model. Section 5.1 describes a Bayesian inference approach using the data augmentation method for general models, including spatial models. Since the data augmentation method is particularly effective for time series models, its application is detailed in Section 5.2.

This study investigated the inverse estimation performance for spatio-temporal data through simulations. Specifically, using the case of the CSNS distribution as an example, we examined whether the estimation performance improves by properly modeling heterogeneity in higher-order moments. To this end, four types of heterogeneity were investigated. The models are defined as follows: the model with $\theta_{tk}=0$ for $t=1,\ldots,T,k=1,\ldots,K$ is referred to as the zero theta (ZT) model, corresponding to the Gaussian regression model; the model with $\theta_{tk}=\theta$ for $t=1,\ldots,T,k=1,\ldots,K$ is referred to as the fixed theta (FT) model; the model with $\theta_{t1}=\cdots=\theta_{tK}=\theta_{t}$ for $t=1,\ldots,T$ is referred to as the temporal varying theta (TT) model; and the model with $\theta_{1k}=\cdots=\theta_{Tk}$ for $k=1,\ldots,K$ is referred to as the spatial varying theta (ST) model.

The observation model is as follows:

where $\bm{\xi}_{\tau}$ follows a CSNS distribution, depending on the $\theta_{tk}$ parameter, using the matrix square root. The covariance matrix $\bm{\Sigma}$, which models spatial correlation, is assumed to follow Equation (3). The time length $T$ for parameter estimation was set to $50$, whereas the time length for evaluating prediction accuracy, $T_{future}$, was set to $2$. The location is represented by a $5\times 5$ grid. The adjacency matrix is defined by assigning a value of $1$ to the neighboring pairs that share a border with separating grids, and $0$ to all others. The features $\bm{x}_{\tau}$ is generated as a random sample from the standard normal distribution.

The true values for the simulation model for each case are as follows: in Case 1, we used the ZT model; in Case 2, the FT model with $\theta=2.5$; in Case 3, the TT model with $\theta_{\tau}\sim N(2.5,3^{2}),t=1,\ldots,T$, where $\theta_{\tau}$ were generated as random variables; and in Case 4, the ST model with $\theta_{k}\sim N(2.5,3^{2}),k=1,\ldots,K$, where $\theta_{k}$ is generated as random variables. For the remaining parameters, the coefficient $\beta_{1}$ for randomly generated features was set to $1.0$, intercept $\beta_{2}$ to $0.5$, standard deviation of the observation noise $\sigma$ to $0.1$, coefficient of the temporal AR model $\rho$ to $0.5$, parameter $\omega$, which determines the variability of $\bm{\theta}$, to $1.0$, and parameter $\eta$, which controls the strength of spatial correlation, to $0.5$.

The prior distributions for the estimation model were shared across all four models, following Lee et al. (2018). We employ a hierarchical distribution as the prior distributions for $\bm{\theta}_{1:T}$. The priors are as follows:

where $\mathcal{U}(0,1)$ denotes the uniform distribution from $0$ to $1$.

Using the sampling method described in Section 4.3, we obtained 1,000 samples from the posterior distribution for each simulation to evaluate estimation performance. During the sampling process, more than 50,000 samples were discarded corresponding to the burn-in period based on the results of trace plots, and thinning was applied to ensure that the autocorrelation was below 0.3. Figure 3 presents an example of a trace plot for several parameters of the FT model in Case 1, after discarding burn-in samples and applying thinning. The plot confirms that convergence was sufficient.

Figure 4 presents the root mean squared errors (RMSE) of each parameter for the four cases and four estimation models that were simulated 50 times. Box plots are used to illustrate each case. In the ZT model, $\bm{\theta}$ was fixed to $\bm{0}_{TK}$ and not estimated; therefore, the RMSE plot for $\bm{\theta}$ is left blank for the ZT model.

First, for Case 1, assuming normally distributed latent variables, the estimation performance is nearly identical across all four models. This indicates that even when the true model follows a normal distribution, accounting for higher-order moments does not degrade estimation performance. For Case 2, the estimation performance of the ZT model is generally poor; however no significant difference in performance is observed for the other three models. This suggests that when $\theta$ is fixed to a single value, consideration of heterogeneity in higher-order moments does not impair estimation performance. Finally, in Cases 3 and 4, the true models (TT and ST) demonstrate improved estimation performance for most parameters. Notably, the improvements are significant for the standard deviation of observational noise ($\sigma$), time-series coefficient ($\rho$), and variability of random effects ($\omega$). These results suggest that modeling heterogeneity in higher-order moments enhances parameter estimation performance. For $\theta$, the FT model performed well in Cases 1 and 2, whereas in Case 3, the variation is slightly different with no notable differences among the three models. In Case 4, despite increased variation, the ST model improved estimation performance.

Similarly, Figure 5 displays predictive accuracy, which was evaluated using log marginal predictive likelihood (LMPL), LMPL in the future data (FLMPL), and RMSE of $y$ in the future data (FRMSE). The calculation method for LMPL followed that used in the CarBayesST package (Lee et al., 2018). FLMPL is an estimator of LMPL for the future data $y_{T+1:T+T_{future}}$ conditioned on the learning data $y_{1:T}$, specifically $\log p(y_{T:T+T_{future}}|y_{1:T})$. FLMPL is calculated from samples of the posterior distribution as follows:

where $\bm{s}^{(i)},\bm{\theta}^{(i)},\bm{\beta}^{(i)},\sigma^{(i)},\rho^{(i)},\omega^{(i)},\eta^{(i)}$ is a sample of parameters from the posterior distribution conditioned on $\bm{y}_{1:T}$, and $S$ is the number of samples. Similarly, FRMSE is the RMSE of $\bm{y}_{T+1:T+T_{future}}$ and is numerically calculated as follows:

where $\hat{\bm{y}}_{T+1:T+T_{future}}^{(i)}$ is a sample drawn from $p(\bm{y}_{T+1:T+T_{future}}|\bm{y}_{1:T})$, and $||\cdot||$ is the Euclidean norm.

For LMPL, in Case 1, there is no significant difference among the four models, whereas in Cases 2, 3, and 4, the true model shows higher values. This indicates that selecting the true model improves the fit to the observed data. In contrast, for FLMPL, the ZT model performs poorly across all cases, with the TT model consistently performing slightly better than the FT and ST models. In this scenario, the standard deviation of the observational noise is 0.1, whereas the FRMSE ranges from approximately 1.1 to 1.5. This can be attributed to the assumption that $\theta$ varies over time, causing the model to exhibit heavier tails, which influences prediction. Finally, for FRMSE, there is little variation across the cases. This can be interpreted as the CSNS distribution maintaining the mean and variance, resulting in minimal differences in squared errors of the observations.

This study focuses on identifying the production functions for 48 states in the United States, using annual data from 1970 to 1986, following the model outlined in Baltagi (2021). A production function represents how capital and labor affect output. We assumed a Cobb-Douglas type production function, where output is expressed as the product of capital and labor raised to respective exponents. By taking the logarithm of the Cobb-Douglas production function and rewriting it in a form that allows for coefficient estimation from observed data, it can be expressd as follows:

Here, $y$ represents the logarithm of the gross state product, with the features including private capital stock ($\mathrm{KP}$), highways and streets ($\mathrm{KH}$), water and sewer facilities ($\mathrm{KW}$), other public buildings and structures ($\mathrm{KO}$), nonagricultural payroll employment ($\mathrm{L}$), and the state unemployment rate ($\mathrm{Unemp}$). All features, except the unemployment rate, are expressed in natural logarithms. The data spanning from 1970 to 1984 were utilized as the training dataset, whereas the data from 1985 to 1986 were used to validate the predictive accuracy of the models. Figure 6 displays the value of $\bm{y}$, the logarithm of the gross state product for each state. Data from the first (1970), middle (1977), and final year (1984) are plotted.

For the models, we compared the estimation results of ZT, FT, TT, and ST, using the same approach and priors as that of the simulation study. We set sufficiently long burn-in and thinning periods and obtained $4000$ samples from the posterior distribution to evaluate the estimation results.

Table 3 displays the median and 90% credible interval of samples drawn from the posterior distributions of parameters for each model. The medians of the posterior distributions of the all parameters exhibit similar patterns, although a slight variation is observed in the estimates of the regression coefficients. This variation can be attributed to the proposed model, which preserves both the mean and variance.

Figure 7 displays the median of the posterior distribution of $\bm{s}_{\tau}$ for each state in 1970, 1977, and 1984 under each model. The figure reveals that, while the trends in the estimates of $\bm{s}_{\tau}$ are similar across the four models, the actual values differ significantly. For example, the posterior median of $s$ for Wyoming in 1970 under the ZT model is approximately $0.15$, whereas it is $0.28$ under the TT model. Certain states appear in much brighter colors under the CSNS model. For instance, Wyoming around 1970 is depicted in a brighter color, indicating higher productivity than expected based on its capital and infrastructure. In 1970, Wyoming was the second least populous state (United States Census Bureau, 1973), exhibiting low feature values, with $L$ being the smallest. However, owing to its abundant resources such as coal, crude oil, and natural gas (U.S. Energy Information Administration, 2023), Wyoming demonstrated a total production that surpassed predictions based on its features. The states appearing in brighter colors under the CSNS model suggest that this model assigns greater probability to positive values than the normal model, because $\theta_{i}$ is positive. Additionally, in the ST model, differences in shading intensity are more distinct separated than in the FT and TT models.

Figure 8 depicts the skewness of the posterior distributions in each model. In this figure, in the ST model for 1970, the differences in skewness are more distinct than in the other models. This is attributed to the ST model’s ability to capture the spatial heterogeneity of higher-order moments. In 1977 and 1984, these differences in skewness are less pronounced, presumably because the variance increases over time, making it more difficult for the skewness to become substantially larger.

Finally, Table 4 presents the evaluation of predictive accuracy. Compared with the ZT model, the FT, TT and ST model show improvements in all evaluated metrics. Even though $\bm{y}$ has undergone a logarithmic transformation, considering the heterogeneity of higher-order moments leads to enhanced estimation performance. In particular, the TT and ST models outperform the FT model, suggesting that incorporating both temporal and spatial heterogeneity in higher-order moments improves model fit and predictive accuracy.

In this study, we propose a method that enables the treatment of the heterogeneity of the higher order moments while preserving the mean and variance, crucial for maintaining parameter interpretability (see Section 1). After demonstrating closure under linear transformations of the proposed model, we analytically derived its skewness and kurtosis. In addition, applying it to spatio-temporal data, we examined how the properties of the model change depending on the chosen matrix decomposition method. In particular, using the matrix square root preserves symmetry with respect to ordering, making it suitable for data with inherent symmetry. For time-series data, we showed that the Cholesky decomposition corresponds to the AR model representation. We then proposed an efficient Bayesian inference method employing data augmentation and the FFBS algorithm.

In the simulation study, we compared the estimation performance of the four models that incorporated heterogeneity in higher-order moments. The results demonstrated that appropriately accounting for higher-order moment heterogeneity improves parameter estimation performance, model fit, and predictive accuracy.

For a real-data application, we applied the CSNS model to identify production functions across the United States. Although parameters other than $\bm{s}$ generally exhibited similar trends, the posterior distribution of $\bm{s}$ revealed that skewness varies substantially by location, particularly when spatial heterogeneity in higher-order moments is considered. Furthermore, both the TT and ST models showed improvements in terms of model fit and predictive performance metrics, indicating that even when the response variable is log-transformed, accounting for higher-order moment heterogeneity enhances estimation performance.

An open question for future research is the development of a model that simultaneously preserves the mean and variance, maintains spatial conditional independence, and accommodates heterogeneity in higher-order moments. Furthermore, as the proposed model encompasses a broad class of models, a potential extension may explore more complex structures where each dimension follows a different distribution, which remains unaddressed in this study.

The Julia implementation of CSNS model is available on GitHub:
(https://github.com/Kuno3/CsnSubclass).

This work was supported by JSPS KAKENHI Grant Number 24K00175.