A flexible model for spatial volatility with an application to the Chicago housing market

ABSTRACT Existing volatility models normally emphasize the behaviour of prices in a temporal sense and comparatively few studies have explicitly analysed the spatial variation of volatility. This paper proposes a flexible spatial volatility model for squared returns using a Box–Cox transformation that includes the linear and log-linear forms as special cases, thus providing a unified framework for simultaneously testing space-varying volatility and its functional form. The use of the model is illustrated by a substantive application to housing price data in the US city of Chicago. The estimation results suggest that housing returns in Chicago show that the volatility exhibits strong spatial dependence and the log-linear functional form is appropriate. In the final log-linear model, a new practical indicator, called neighbourhood elasticity, is proposed that determines how volatility in one neighbourhood is linked to that in surrounding neighbourhoods.


INTRODUCTION
The volatility of housing prices has important implications for household behaviour and welfare as well as for the aggregate US economy. At the household level, it can be easily argued that housing is the most important asset for many households. It is usually both the largest asset they own and the most readily available source of collateral against which they can borrow. 1 Higher housing price volatility thus has the potential to pose substantial risk to household welfare. For example, it could distort a household's housing choices, lead to a higher likelihood of mortgage foreclosure, and also affect home building and intergenerational equity (Miller & Peng, 2006; Oxley & Haffner, 2010; Stephens, 2011. From a macroeconomic perspective, the impact of housing price volatility is similarly damaging as the housing sector is vital to the national economy. Recent experience has made painfully clear the importance of the housing market in the United States. A catastrophic and systematic collapse of the US housing market triggered an economic recession, the so-called Great Recession, which rippled throughout the global economy.
Research into housing price volatility has received increased attention in the housing literature in recent years. Much of this research has focused primarily on investigating whether housing price volatility is time-varying, that is, housing prices exhibit volatility clustering or autoregressive conditional heteroskedasticity (ARCH) effects (e.g., Barros et al., 2015; Crawford & Fratantoni, 2003; Dolde & Tirtiroglu, 1997; Miles, 2008; Miller & Peng, 2006. While it is widely recognized that conditional heteroskedasticity is pervasive in studies of housing price volatility, there has been little research addressing such heteroskedasticity in the context of spatial volatility.
In order to see the spatial variation of volatility in the local housing market in the United States, consider the annual return on housing in 2017 at the census tract level in Chicago, Illinois. Construction of the return data is described in section 5. Looking at the representative spatial plot for the returns and squared returns as a volatility proxy in Figure 1, it is evident that while returns appear to be randomly distributed over space, squared returns are spatially correlated, with distinct clusters of high and low volatility readily identifiable, for example, high volatility clusters are detected in Chicago's west-side and south-side neighbourhoods. 2 This clustering pattern may occur in housing markets in large cities, as the nature of change in the US housing market has opened significant gaps between regions in recent decades. As suggested by one of the referees, I supplemented Figure 1 with the box plots in Figure 2.
In addition to the spatial quantile maps and box plots, the local indicators of spatial association (LISA) are also used to explore the spatial clusters and presented in Figure 3. The analysis indicates that low volatility is geographically clustered in the north-east side while high volatility is clustered in the west and south sides of the city.
Motivated by the visual implications presented above, this paper develops a spatial ARCHtype model of squared returns to examine housing price volatility in the context of a spatial regression framework. 3 Specifically, the proposed model incorporates space-varying volatility (or spatial volatility clustering), which is the spatial equivalent of time-varying volatility in a time-series context.
When analysing economic data that take only positive values, the use of transformations is very common and may be helpful when the usual assumptions are not satisfied in one's original model. In the finance literature, the logarithmic transformation of realized volatility (the sum of squared intraday returns) is often used in empirical applications owing to its superior finite sample properties. Recently, however, a more flexible power transformation suggested by Box and Cox (1964) has been considered in this context (Gonçalves & Meddahi, 2011; Nugroho & Morimoto, 2016; Taylor, 2017; Weigand, 2014; Zheng & Song, 2014. The Box-Cox transformation indexed by the transformation parameter defines a general class of functional forms that includes the linear and log-linear forms as special cases. This feature allows the data added flexibility in model specification and provides a unified structure for statistically distinguishing between alternative functional specifications. In light of the foregoing, this paper generalizes the spatial ARCH model to a more general version with the Box-Cox transformation to derive the most appropriate functional form of spatial volatility, where the linear and log-linear models are special cases (hereafter referred to as the BC-SARCH model). The advantages of the BC-SARCH model are that both functional form and spatial dependence can be considered simultaneously, and that the validity of linear or log-linear specifications can be addressed. The proposed model would also be useful in another respect, in that the transformation might affect the distributional shape of the variable favourably. It is well known that the unconditional distribution of data for which typical linear models are used is frequently skewed and leptokurtic. Given that appropriate transformation may induce  symmetry to the distribution, it seems probable that the proposed model would yield, in addition to capturing the presence of spatial dependence, approximate symmetric and mesokurtic properties (i.e., normality) for the unconditional distribution (Sarkar, 2000).
To estimate the BC-SARCH model, a maximum likelihood estimation (MLE) procedure is applied to both simulated and empirical data. The finite sample properties of the estimator is examined in a Monte Carlo experiment and the result suggests that the MLEs perform well in estimating parameters of interest. This study also carries out a detailed empirical analysis using housing sales price data within the city of Chicago in a period running from 2009 to 2018. In this empirical application of the BC-SARCH model, the results suggest that volatility exhibits substantial spatial dependence and its functional form is close to the log-linear model. Judging by the associated model diagnostics and specification tests, the log-linear model is found to be superior to the linear model. Furthermore, the resulting spatial dependence parameter in the log-linear model serves as neighbourhood elasticity, which measures how volatility in one neighbourhood is linked to that in surrounding neighbourhoods. This analysis is further extended to the spatial panel data models to test temporal heterogeneity in neighbourhood elasticity, that is, whether neighbourhood elasticity stays the same over time. Applying the method introduced by Xu and Yang (2020), this paper shows that, after controlling for both spatial and temporal heterogeneity in intercepts, neighbourhood elasticity is homogeneous and there exists no structural change.
This paper contributes to the existing literature methodologically and empirically. The methodological contribution consists of a unified structure for statistically testing both spatial dependence and functional form simultaneously. The empirical contribution relates to the use of the most recent multiple listing service (MLS) price data, allowing extension of spatial volatility analysis at a finer spatial scale (the census-tract scale). As noted above, while the spatial dependence of housing price volatility has been detected at large spatial scales (the state or metropolitan statistical area (MSA) scale), 4 a rigorous empirical examination of the linkage at a much smaller spatial scale has not been carried out because of the lack of publicly available data. 5 Finally, a new practical indicator, neighbourhood elasticity, is proposed that can serve as a tool for policymakers to assist them to mitigate volatility transmission and the risk of contagion in the housing market.
The rest of this paper is organized as follows. In Section 2 I review the related literature on spatial volatility models. In Section 3 I discuss the detailed description of the BC-SARCH model under consideration. The Monte Carlo simulation experiment is conducted in Section 4 to evaluate the finite sample performance of the MLEs for the proposed model. In Section 5 I describe the data used in this paper and associated estimation results. Section 6 concludes with a discussion of the contributions and future developments of this paper. The robustness of the findings in alternative modeling choices, as well as additional figures and tables can be found in the Appendix in the supplemental data online.

RELATED LITERATURE
An increasing number of studies have attempted to model the housing price volatility of individual housing markets by employing the class of generalized autoregressive conditional heteroskedasticity (GARCH) models. Earlier studies have employed different univariate GARCH-type models, such as GARCH and GARCH-in-Mean (GARCH-M) (e.g., Dolde & Tirtiroglu, 1997; Hossain & Latif, 2009; Stevenson et al., 2007, exponential GARCH (EGARCH) and exponential GARCH-in-Mean (EGARCH-M) (e.g., Lee, 2009; Lin & Fuerst, 2014; Willcocks, 2010, component GARCH (CGARCH) and component GARCH-in-Mean (CGARCHM) (e.g., Karoglou et al., 2013; Lee & Reed, 2014; Miles, 2011. Nevertheless, there remains a lack of statistical models accounting for the spatial conditional heteroskedasticity in a spatial econometric context, despite the growing body of evidence that neighbourhood plays an important role in volatility dynamics as mentioned in the introduction. Subsequent to the pioneering work of Bera and Simlai (2005), who propose a special type of spatial ARCH (SARCH) model by employing the information matrix (IM) test statistic in a simple spatial autoregressive (SAR) model, several ARCH or GARCH-type models have been proposed to incorporate both temporal and spatial effects in the dynamics of volatility. For example, Borovkova and Lopuhaa (2012) and Caporin and Paruolo (2006) introduce a temporal GARCH model, which includes temporal lags influenced by neighbouring observations, while Otto et al. (2018) introduce an exponential SARCH model and illustrate the use of model as a residual process for the spatial modeling of lung cancer mortality in US counties. Very recently, Sato and Matsuda (2021) suggest a further extension of SARCH models with applications to land prices in the Tokyo area of Japan.
In addition to these time-series based models, Le Gallo et al. (2020) develop tests to detect spatial group-wise heteroskedasticity (SGWH) based on the scan methodology in the SAR and spatial error model residuals and use them to explain clustering in the Madrid housing prices for the attics submarket. 6 Several other recent papers have discussed moment methods of spatial econometric models with heteroskedasticity (e.g., Breitung & Wigger, 2018; Taşpınar et al., 2019.
Interestingly, in spite of this promising start, the use of a spatial econometric approach to detect the presence of conditional spatial volatility has remained largely unexplored in the literature. Furthermore, it is not yet established which functional form is the most appropriate for evaluating space-varying volatility in the housing market.

THE MODEL AND ECONOMETRIC METHODOLOGY
The most widely used model for describing nonlinear dependencies of returns in a time-series context is the ARCH class of models, where past squared returns are used to predict subsequent squared returns. This approach potentially suggests that, in a spatial context, squared returns can be described by the squares of neighbouring observations. The simplest specification for approximating this expectation relates the squared returns in one location to the weighted average of the neighbouring squared observations. To specify the functional relationship, the following BC-SARCH model is proposed, which is expressed as: where y 2 i is squared returns at location i, y 2 j is the values of the variable at surrounding locations j (j = i), w ij is the spatial weights, which are non-zero when i and j are neighbours, and zero otherwise, a 1 is the spatial dependence parameter, l is the transformation parameter, and e i IIDN (0, s 2 ), i = 1, 2, . . . , n. It is important to note that no exogenous variables are introduced in this model as the data on typical hedonic pricing variables (e.g., structural or neighbourhood attributes of housing) were not sufficient to attain statistical significance in my previous analysis. 7 Equation (1) can be also expressed using the matrix form by letting y 2 = (y 2 1 , y 2 2 , . . . , y 2 n ) ′ be the data on which the Box-Cox transformation is applied: where y 2(l) is the Box-Cox transformation of y 2 , W is a row-normalized spatial weight matrix such that n j=i w ij = 1 for all i, with zeros on the main diagonal, the off-diagonal elements take values between 0 and 1, and 1 is an n-dimensional vector of ones.
The Box-Cox model in equation (1) clearly reduces to the linear model when l = 1, and to the log-linear model when l 0, such that: and thus tests of l aid in selecting the proper functional form. Note that the parameter estimate of a 1 in equation (4) amounts to an elasticity, but one that is conditional on neighbourhood information N. To see this, let the conditional expectation of equation (4) be given by w ij logy 2 j , the neighbourhood average. This in turn leads to: which measures the sensitivity of the expected volatility at location i to an increase in volatility of the neighbourhood, and it will be convenient to call it the neighbourhood elasticity.
The process of choosing the appropriate functional form for the BC-SARCH models involves maximizing a log-likelihood function. The likelihood function is expressed by: where u is a vector of parameters, e = (I − a 1 W )y 2(l) − a 0 1 with 1 denoting the n × n identity matrix, and J is the Jacobian of the transformation from e to y 2 , or: Consequently, the log-likelihood function is given by: where the v i 's are the eigenvalues of W ; see Anselin (1988). By assumption, the matrix I − a 1 W is non-singular. The last term in equation (8) vanishes in the case of the linear model. Finally, the log-likelihood function is maximized with respect to the parameter vector u = {a 0 , a 1 , l, s 2 } to obtain the maximum likelihood estimators (MLEs). The explicit forms of the first and second derivative of the log-likelihood are given in Section A in the Appendix in the supplemental data online. It is important to note that the asymptotic distribution of the MLEs is not formally established in the paper but is likely to hold under similar sets of assumptions developed by Lee (2004), who presents a comprehensive investigation of the asymptotic properties of the MLEs widely used in the literature to estimate spatial models.

Experimental design
In this section I present a Monte Carlo experiment to evaluate the empirical evidence pertaining to the spatial regression using the BC-SARCH specification in a simulated setting. It has been argued in the literature that a spatial regression with heteroskedastic innovation can lead to significantly biased ML estimators (Kelejian & Prucha, 2010). A recent study by Piras and Prucha (2014) points out that the finite sample properties are crucial for understanding bias for estimators of various spatial models. Following the lead provided by the existing literature, I investigate the finite sample performance of the BC-SARCH estimators of interest.
In all the experiments, the data-generating process is assumed to be of the form specified in equation (1). The spatial weight matrix W is generated according to both queen and rook criteria on regular m × m grids, leading to a sample size of n = m 2 . While the queen specification states that two polygons are neighbours if they share a border or a vertex, the rook specification states that two polygons are neighbours if they share a border. The contiguity matrix is typically used as the spatial weight matrix for data represented by areal units (polygons) that vary in size. Both weight matrices are row-normalized so that each row sums to unity. This study considers all combinations of the true parameters a 1 [ {0.3, 0.5, 0.8} and l [ {0.1, 0.3, 0.5, 0.8, 1.0}, and the values of a 0 and s 2 are set equal to 1.0 and 0.1, respectively. Note that the parameter space of a 1 is specified such that I − a 1 W is non-singular, which requires that where v min and v max , respectively, are the smallest (i.e., the most negative) and the largest real eigenvalues of W . For each combination of a 1 and l, this procedure is repeated 1000 times using samples of size n = 100, 400 and 900. Following Kelejian and Prucha (1999), the quality of an estimate is evaluated in terms of its bias and the root mean squared error (RMSE). The measure of bias is defined as the average difference between the estimator and the true value, while the RMSE is defined as the square root of the sum of the variance and the squared bias of the estimator.

Simulation results
The simulation results using the queen contiguity are reported in Table 1. The first two columns are true values of the parameters of interest, a 1 and l, and the remaining columns present the estimated bias and the RMSE of the two parameters for various combinations of sample size n and parameter values a 1 and l.
The results suggest that, on average, the bias of a 1 does not follow any particular pattern with respect to the parameter set of high or low values of a 1 and l. Nonetheless, for a large true value of a 1 , a 1 = 0.8 here, the estimator of a 1 displays small bias, and the average RMSE is the lowest across all values of l. Especially for a 1 = 0.8 and l = 0.1, the estimator of a 1 exhibits the lowest RMSE. Comparing the results associated with each value of l across all sample sizes, it is found that increasing the sample size improves the quality of the estimates, that is, the bias and RMSE become uniformly smaller with larger samples.
Another important parameter of interest is l, which is the transformation parameter in the BC-SARCH model. It is shown that the estimator performs reasonably well in estimating l in terms of bias and RMSE. In general, for small true values of l, the estimator tends to produce better estimates with small bias and RMSE. In particular, the estimator with l = 0.1 has the smallest RMSE across all values of a 1 and sample sizes, while the largest bias corresponds to a situation where the model is a linear form (l = 1). The RMSE is shown to decrease as the sample size increases. The simulation results using the rook contiguity matrix are in line with the findings from the results based on the queen contiguity, suggesting that the performances Table 1. Bias and root mean square error (RMSE) for estimators of α 1 and λ based on queen contiguity. n = 100 n = 400 n = 900 of the estimators are robust to spatial weight matrix specifications (results available from the author upon request).

Data
I base the empirical analysis on residential sales data obtained through the multiple listing service (MLS) used by Illinois REALTORS ® . 8 The data include records of all residential property sales in the city of Chicago and cover the 10-year period running from 2009 to 2018. These properties are comprised of single-family properties and condominiums. To remove outliers, transactions that take place at extreme prices, below the first or above the 99th percentile of the distribution of raw prices, are excluded. I chose census tracts as the spatial units for the analysis, owing mainly to their homogeneity with respect to socio-economic and demographic characteristics. 9 Furthermore, most census data are reported at this level of geography and the boundaries of the census tracts follow permanent and easily recognizable physical features.
There are 801 census tracts in Chicago and, for ease of interpretation, 77 community areas are used as references for the neighbourhoods. 10 A detailed description of the areas can be found in Table D1 in the Appendix in the supplemental data online. Figure 4 shows a map of Chicago with the 801 census tracts outlined in light blue and the 77 community areas outlined in dashed black. It should be noted, however, that several tracts in which there are no observed transactions are eliminated from the analysis. This led to reduced and varying sample sizes across years about 759-776 observations (see Table D2 in the Appendix online).  Table D1 in the Appendix in the supplemental data online for details). Source: US Census Bureau, OpenStreetMap.
All transaction records within the study area are geocoded using ArcGIS 10.6 and aggregated to the census-tract level to calculate an annual median price for each tract. The rationale for using median home prices in this study is its long history and extensive coverage of census tracts. The calculated median prices are then used to construct the key variable of interest in this study, housing returns. The return series is obtained from the first difference of log annual median house prices, that is, y i,t = log(P i,t /P i,t−1 ), where at census tract (location) i, P i,t and P i,t−1 are the median prices at time t and t − 1, respectively. 11 It should be noted that following Willcocks (2010), the returns are not smoothed or adjusted for inflation as adjustment could hide the impact of volatility changes between adjacent time periods.
To investigate the spatial dependence structure in volatility, the series of squared returns y 2 i and log-squared returns log y 2 i are considered as proxies of the volatility measures, which are the linear and log-linear cases nested within the BC-SARCH model. The logarithm of squared returns is problematic when the returns are zero or very small in magnitude. To avoid the problem of taking the logarithm of zero, Fuller (2009) proposes a slight modification of the squared returns, which is given by: where s 2 is the sample variance of returns and g is a small constant, with g set equal to 0.02, following Fuller. This modified version will, for convenience, be hereinafter referred to as the logsquared returns. Summary statistics for returns and two volatility series, squared returns and log-squared returns, are presented in Table D2 in the Appendix in the supplemental data online, including the Jarque-Bera (J-B) test statistics for normality for each year during the study period running from 2009 to 2018.

Preliminary tests
To understand whether the three return series reflect spatial dependence, a global Moran's I test for spatial patterns is undertaken for cross-sectional data from 2009 to 2018 separately. The Moran's I is a widely used test for the presence of spatial dependence, and in this study the queen and rook contiguity weight matrices are considered to check the robustness of the test results. Table 2 provides the test results based on the queen contiguity matrix. Note that results  based on the rook contiguity matrix are very similar as those based on the queen contiguity matrix and available from the author upon request. The results are in line with expectations presented visually in the introduction and are not sensitive to the choice of the spatial weight matrix (results are available from the author upon request). Both measures of volatility exhibit a high degree of spatial dependence, while there is no or only weak statistical evidence for spatial dependence in the returns, with the exception for 2009, when the effects of the financial crisis led to significant negative returns for most tracts, which can be identified in panel (a) of Figure 5. 12 On the whole, the above results appear to be analogous to the well-known stylized fact in the time-series literature that, in contrast to the lack of serial dependence in returns, the autocorrelation for the squared returns is always positive and significant. Note that, in the case of log-squared returns, the reported test statistics are generally larger, and p-values are lower. The presence of spatial dependence in the return series indicates that the housing market is spatially inefficient, and therefore the housing returns could contain useful information for predicting the returns in neighbouring locations. This result is not surprising because the unique features of a housing market, such as transaction costs, infrequent transactions, and a high degree of heterogeneity, may limit arbitrage opportunities, leading to pricing inefficiencies (Case & Shiller, 1988).

Cross-sectional estimation results
The BC-SARCH model specified in equation (1) is estimated separately for each of the ten years beginning in 2009 based on the queen contiguity matrix, and the cross-sectional results are   Table 3, where the parameters are estimated by the maximum likelihood method. After estimating the coefficients, their corresponding standard errors are obtained from the diagonal elements of the inverted information matrix and used to perform a t-test for the null hypotheses that each of a 0 , a 1 , s and l is equal to zero.
The most noticeable result reported in Table 3 is that, as expected, there is strong and significant spatial dependence in Box-Cox transformed squared returns in all years, with coefficient values ranging from 0.293 for 2013 to 0.613 for 2009. This result provides strong evidence for the presence of spatial clustering in volatility discussed earlier in connection with Figure 1 and Table 2. That is, locations with high volatility tend to have neighbours with high volatility and those with low volatility tend to have neighbours with low volatility.
Of particular interest is the transformation parameter l. In every year, the estimated coefficients appear to be very close to zero, ranging from 0.094 for 2010 to 0.132 for 2013. The t-test indicates that the null hypothesis that l = 0 is not rejected and thus l is not significantly different from zero, pointing to a log-linear form as the best fit. If the test is based on the null hypothesis that l = 1, the null hypothesis is clearly rejected for all years (results available from the author upon request). Therefore, the estimation results of the BC-SARCH model applied to the empirical data indicate that there is strong evidence against a linear model and in favour of a log-linear model with l 0. Finally, at the bottoms of each panel in Table 3, several model specification tests are performed to verify the validity of the estimated BC-SARCH model using Moran's test, Lagrange multiplier (LM) test and the J-B normality test for the residuals. The Moran's I test as described by Kelejian and Prucha (2001) 13 fails to reject the null hypothesis of the absence of spatial dependence. Additionally, the LM test does not reject the same null hypothesis. Therefore, these results suggest that the spatial dependence has been entirely captured by the model and is not left in the residuals. 14 With regard to normality, J-B statistic confirms the normality behaviour of the estimated residuals. The same conclusion can be drawn when both estimation and specification tests are replicated using alternative spatial matrices, which further strengthens the usefulness of the BC-SARCH model. The robustness analysis can be found in Section C in the Appendix in the supplemental data online.
All in all, the above results support that the presence of conditional volatility over space is non-trivial and the proposed BC-SARCH model is a powerful tool for modeling such spatially dependent heterogeneity. The results further reveal that a model of the log-linear form may provide an adequate representation of spatial volatility. From an economic perspective, the most well-known benefits of using a log-linear model are that it cannot only approximate a normal distribution and reduce skewness but also make it easier to interpret the empirical results. While a Box-Cox transformation severely complicates the interpretation of the coefficients, the coefficients of a log-linear model lend themselves to straightforward interpretations such as elasticity, a measurement of the responsiveness of one variable to a change in another variable in percentage terms. Therefore, as discussed in Section 3, a log-linear model allows a 1 to be interpreted as neighbourhood elasticity, which measures the percentage change in volatility of neighbouring locations with respect to a one percent change in volatility in one location. Table 4 reports the estimation results of the final log-linear model in equation (4), which are obtained using the ML method. The results reveal that, despite the simplicity and convenience of the model, the main findings remain unaltered. The estimated coefficients of the spatial dependence parameter a 1 , that is, neighbourhood elasticity, are still positive and highly significant at the 1% significance level across all years with almost the same magnitudes as those obtained with the BC-SARCH model, suggesting an inherent volatility dependence among cross-sectional units. Values of the coefficients range from 0.296 for 2013 to 0.623 for 2019 with a cross-sectional average of 0.435 for the case of queen contiguity, which suggests that a  (Anselin, 1988). Diagnostic tests, including Moran's I, Lagrange multiplier (LM) for spatial lag model and Jarque-Bera (J-B), are implemented to examine the spatial dependence and normality of the residuals. * * * p < 0.01; * * p < 0.05; * p < 0.1. Note: Standard errors are shown in parentheses. Pseudo-R 2 is computed as the squared correlation between the observed and predicted values of the dependent variable (Anselin, 1988). Diagnostic tests, including Moran's I, Lagrange multiplier (LM) for spatial lag model and Jarque-Bera (J-B), are implemented to examine the spatial dependence and normality of the residuals. * * * p < 0.01; * * p < 0.05; * p < 0.1.
1% increase in volatility leads to an approximately 0.4% increase in volatility in neighbouring locations.
Conducting the same diagnostic tests as I applied to the BC-SARCH model for the log-linear model reveals a similar regime. Moran's I test confirms that no evidence of remaining residual spatial dependence is found. Although the J-B test rejects normality, the small values of the statistic imply that there is no severe departure from normality.
To further examine how well the log-linear model identifies the spatially dependent pattern of volatility, the estimated squared returns obtained from the model based on 2017 and the queen contiguity are plotted over the census tracts on the three-dimensional (3D) maps with two views presented in Figure 6. As can be seen, the spatial dependence pattern is well reproduced by the model. While the lowest estimates of squared returns are distributed over the north-east side past downtown Chicago, the estimated squared returns show clear spikes on the west and south sides of the city, with distinct clusters. All in all, the simplified log-linear model appears to be an appropriate tool with which to describe the spatial dependence of volatility in the Chicago housing market. 15

Pooling the time-series and cross-section data
Although the cross-sectional results show consistent evidence of the presence of spatially dependent volatility under multiple model specifications, the estimates point to the need to avoid unreliable conclusions by appropriately accounting for the time-varying heterogeneity that arises from unobserved factors to avoid unreliable conclusions. The rationale is provided in Figure 7, which plots the estimated time-varying spatial coefficients of α 1 together with the 95% pointwise confidence intervals obtained from the BC-SARCH and the final log-linear models over the 2009-18 period, both based on the queen contiguity matrix. Figure 7 shows that, over the entire observed time period, the coefficient of a 1 was substantial and almost identical across the two models, indicating that neighbourhood plays an important role in explaining volatility in the Chicago housing market. The magnitude of this effect, however, is not constant over time.
The highest values of the coefficient were observed in 2009 after the outbreak of the financial crisis, after which values started to decrease in 2010 and reached their lowest level in 2013. After 2013, the coefficient exhibits a continuously increasing trend. In summary, the results indicate the varying but moderate persistence of the spatially dependent volatility dynamics over the last decade in the Chicago housing market.
Several studies in the housing literature suggest that the degree to which volatility depends on spatial dependence varies across market conditions. For example, Miao et al. (2011) show that Figure 6. Estimated squared returns based on queen contiguity, 2017. spatial linkages appear more intense during the active phase (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006) than during the tranquil phase . Zhu et al. (2013) further define the period after 2006 as the crisis phase (2007-09) and observe the increased intensity during the crisis period. Moreover, they find that, compared with spatial interactions that occur during tranquil and boom phases, spatial interdependence is much more distinct during the crisis. The present study suggests that such patterns identified at a large spatial scale (the MSA scale) remain similar at a much finer spatial scale (the census-tract scale). While the strongest spatial dependence of volatility between regions is detected during the financial crisis, a sizable increase in the dependence is also observed during the recent upturn that began in 2013 when housing prices started to recover from the declines recorded during the crisis. In this sense, the aspects of time-varying heterogeneity across census tracts have to be considered in the estimation and adding the time-varying spatial effects to the regression model may be useful.
The literature on temporal heterogeneity in spatial panel data models is, however, limited. 16 Only a handful of studies have considered structural breaks in heterogeneous spatial panels. Sengupta (2017) considers hypothesis-testing for a structural change in a spatial panel model without fixed effects, while Li (2018) proposes and studies fixed-effects spatial panel data models with structural change using the quasi-maximum likelihood method. Angulo et al. (2017) focus on the problem of testing for changes in the weight matrix over time.
Some scholars in the finance literature have also attempted to make similar explorations of spatial panel data models with time-varying parameters. For example, Blasques et al. (2016) and Catania and Billé (2017) extend the spatial static panel data model by introducing time-varying spatial dependence, but their models are most suitable for a case with a large time dimension and a small cross-sectional dimension.
Very recently, Xu and Yang (2020) consider fixed-effects spatial panel data models with temporal heterogeneity in regressions and spatial coefficients by focusing on testing problems. They introduce a general method, the adjusted quasi-score (AQS) method, for constructing the specification tests for temporal homogeneity/heterogeneity in regression and spatial coefficients in spatial panel data models, allowing for the existence of spatial and temporal heterogeneity in the intercepts (or fixed effects) of the model. This study presents the first application of their method to empirically test temporal heterogeneity in the spatial dependence of volatility, neighbourhood elasticity in the log-linear model, in the housing market.
Changing notations to adapt to the changing context in which a dependent variable is logsquared returns and no explanatory variable is included as described in equation (4), first consider the following simplest panel spatial lag model with individual specific fixed effects or one-way fixed effects (1FE-SL): where a 1t is neighbourhood elasticity at time t, u i denotes the individual-specific fixed effects or the spatial heterogeneity in the intercept, and z it is an n × 1 vector of i.i.d. disturbances with mean zero and variance s 2 . As the aim of the present study is to test the hypothesis that neighbourhood elasticity stays the same over time, only tests of temporal homogeneity in spatial coefficients, that is, H 0 :a 11 = · · · = a 1T = a 1 , are conducted, allowing for the presence of unobserved cross-sectional heterogeneity in intercepts, that is, the individual fixed effects u i . They further extend their tests to a panel spatial lag model with both individual and timespecific fixed effects or two-way fixed effects (2FE-SL), which can be written in the following form: where v t are the unobserved time-specific effects or the unobserved temporal heterogeneity in the intercept. Based on these models, two types of AQS tests (naive and robust) are proposed and these tests are fully extended to the spatial panel data models with both spatial lag and error (SLE) dependence, where the disturbances are also subject to spatial interactions. However, since they show that the robust tests have much superior finite and large sample properties than the native tests through Monte Carlo simulation, and a model specification search is not of primary interest in the present study, only asymptotically valid and non-normality robust AQS tests and spatial lag models given in equations (10) and (11) are considered.
As the proposed tests are intended for balanced panel data, only census tracts that have the same number of time observations are included to obtain a balanced panel. The purpose of constructing the pseudo-panel data is to model time heterogeneity and make full use of information in the data, improving the accuracy and efficiency of estimation. 17 The resulting pseudo-panel consists of a total of 716 census tracts over 10 time periods between 2009 and 2018. The special weight matrix is based on the queen contiguity matrix and then row-normalized as before.
Returning to Figure 7, it can be suggested that around 2013 there is a shift in the estimated value of neighbourhood elasticity a 1 . Thus, the tests are conducted on the full period  and two subperiods (2009-13 and 2014-18). The 2009-13 period is further broken down into two subperiods, 2009-10 and 2011-13, to examine whether the structure changes during the initial recovery phase after the financial crisis. Table 5 summarizes the values of the test statistics and their p-values for the robust AQS tests for temporal homogeneity based on both the full period and several subperiods, fitted using the two models: a one-way fixed-effects spatial lag model (1FE-SL) and a two-way fixed-effects spatial lag model (2FE-SL).
From the abovementioned results, it can be seen that when the 1FE-SL is applied, that is, temporal heterogeneity in intercepts is not controlled for, tests based on the full data clearly reject the null hypothesis of temporal homogeneity in neighbourhood elasticity. To test whether the null hypothesis holds when 2013 is a change point, the same set of tests is applied to the two subperiods, 2009-13 and 2014-18. The resulting statistics reject the null hypothesis and exhibit no change point. The results for 2009-10 and 2011-13 are also fairly stable, consistently rejecting the null hypothesis. Overall, a rejection of the null hypothesis in any subperiod indicates that neighbourhood elasticity is temporally heterogeneous and it seems to stem from full heterogeneity instead of the existence of change points.
When the 2FE-SL model is applied, however, the null hypothesis of homogeneity of neighbourhood elasticity is not rejected in either full period or any subperiod. This suggests that, when the time fixed effect is taken into account, neighbourhood elasticity becomes homogeneous. It is important to note that when the null hypothesis of homogeneity is accepted in the full period, tests for subperiods are redundant in the sense that acceptance of the null hypothesis for full period also results in accepting the null hypothesis for any subperiod. Nonetheless, the test results for the subperiods lend further support to the homogeneity of neighbourhood elasticity. These findings hold irrespective of the spatial weight matrix used (results available from the author upon request). In summary, the test results show that, after controlling for both spatial and temporal heterogeneity in the intercepts, the homogeneity of neighbourhood elasticity is achieved and there exists no structural change in the term, and moreover a homogeneous spatial panel data model can be used and specifying change points is not necessary.

CONCLUSIONS
This paper studies spatial variation of volatility in the Chicago housing market by proposing a flexible spatial volatility model for squared returns based on a Box-Cox transformation, a technique that has been frequently used as both a flexible functional form and as a decision device to distinguish between alternative model specifications.
The estimation results suggest that housing returns in Chicago show strong spatial dependence in volatility and the commonly used log-linear functional form is appropriate, irrespective of variations in neighbourhood criteria. The appropriateness of the log-linear model is also determined through associated model diagnostics and specification tests. This result has important economic implications for both researchers and housing market practitioners. In the final log-linear model, a new practical indicator is proposed. Neighbourhood elasticity, captured by the parameter a 1 in the model, provides a measure of how volatility in one neighbourhood is linked to that in surrounding neighbourhoods. The average annual elasticity is found to be 0.4, which can be used as a benchmark for comparing housing markets. The findings have important practical implications. The spatial volatility clustering can be used in the forecasting context to obtain appropriate confidence interval and thus helps policymakers develop strategies to mitigate the impacts of volatility transmission and the risk of contagion in the housing market. Finally, to identify whether neighbourhood elasticity remains constant over time, AQS tests for testing the presence of temporal heterogeneity in spatial parameters in spatial panel data models are considered. The test results reveal that neighbourhood elasticity becomes homogeneous after controlling for both spatial and temporal heterogeneity in the intercepts of the model. This paper can be extended by (1) investigating the size and power of a test on the coefficient estimates in the simulation experiments, (2) testing for the existence of spatial non-stationarity in the mean and compute its spatial correlogram for several neighbourhood, and (3) testing for the existence of spatial group-wise heteroskedasticity in the error terms by estimating a hedonic price model. These extensions are important but beyond the scope of the current paper, thus will be addressed in my future studies. Bollerslev et al. (1992), squared returns of not only exchange rate data but also all speculative price series typically exhibit volatility clustering and ARCH-type models are appropriate for volatility estimations. 4. Lee (2009) examines the determinants of housing price volatility for eight capital cities in Australia by using consumer price index, income, unemployment rate and mortgage rate. Miller and Peng (2006) also investigate interactions between housing volatilities and economic conditions in 277 MSAs. They demonstrate that the volatility series is Granger-caused by the home appreciation rate and the growth in gross metropolitan product (GMP). 5. Although Bera and Simlai (2005) use data from 506 census tracts in the Boston metropolitan area, the data are taken from the 1970 census, which does not reflect recent housing cycles. 6. Although Le Gallo et al. (2020) and the present paper have similarities in a sense that both try to detect spatial clusters of observations, they focus on model residuals rather than on the model itself. Therefore, an important difference from Le Gallo et al. is that rather than assuming a certain type of model specification such as SAR and then investigating their residuals, I allow for added flexibility in determining the appropriate degree of non-linearity by developing a new type of model. 7. An editor suggested using the quadratic terms of the explanatory variables to explain the volatility of prices in different locations. However, I found the corresponding parameter estimates to be insignificant. 8. The MLS is a private real estate-listing service where real estate property information is listed and searched by participating members (e.g., real estate agents). 9. Census tracts are small, relatively stable spatial units with population ranging between 2500 and 8000, with an average of approximately 4000, and designed to be homogeneous with respect to population characteristics, economic status and living conditions. 10. The community areas are well-defined, stable geographical regions designated by the Social Science Research Committee at the University of Chicago and officially recognized by the city of Chicago. 11. Hayunga et al. (2019) argue that homeowner characteristics influence their maintenance and home-improvement behaviours, which in turn affect home values and thus the observed pecuniary return. This analysis requires mortgage information data at initial purchase, which are not available from the MLS records; the current study thus treats observed housing returns as actual returns despite the possibility of overstatement. 12. In the finance literature, several papers have documented that an autocorrelation of stock returns exists, particularly during the postcrisis period (Islam et al., 2007; Sarwar & Khan, 2017. 13. Kelejian and Prucha (2001) introduce a central limit theorem (CLT) that can be used to establish the asymptotic distribution of the Moran's I test under certain regularity conditions. See Kelejian and Prucha's equations (2.4) and (4.1) for the test statistic and asymptotic distribution, respectively. 14. The LM test in Table 3 is based on spatial lag model. Other LM test statistics are consistent and available from the author upon request. 15. The fact that the log-linear model is less misspecified than the linear model can also be illustrated by several test statistics based on the information matrix as described in Section B in the Appendix in the supplemental data online. 16. Millo and Piras (2012) develop the splm package in R to estimate the fixed effect spatial panel model that accounts for both spatial lag and spatial error effects, but it does not take account of temporal heterogeneity. 17. Most existing spatial panel estimation methods are designed for balanced panel data, and are not effective for unbalanced panels because of the computational burden associated with inverting a large spatial weighting matrix.