Semi-parameterized Three-dimensional W Matrix

(1)

Semi-parameterized Three-dimensional W Matrix

Chang Tan

∗

S3043800

August 11, 2019

Abstract

The spatial weight matrix W reflecting the geographical, temporal or economic relationships among units is a fundamental part of each spatial econometric model. There are different ways to specify W, but one methodological step forward would be to estimate its elements rather than to specify them in advance. Based on that premise, we propose a new semi-parameterized three-dimensional approach: the strength of distance and time decay effects are parameterized and estimated using non-linear estimation techniques within a SAR model framework, while the degree of similarity with other units is derived from pre-specified in-formation. Furthermore, by utilising the property in real estate models that the matrix W is lower-triangular, we are able to circumvent the perfect solution problem first detailed in Halleck Vega & Elhorst (2015) for spatial econometric models other than the SLX model. Housing transaction data in the Dutch city of Groningen over the period 1993-2017 is used to construct and estimate the proposed semi-parameterized three-dimensional W matrix, while Moran’s I is applied to evaluate its performance.

1 Introduction

Spatial econometric models are applied to study the impact of spatial lags and to measure spatial spillover effects among different units. Both effects are based on specifications of the

(2)

W matrix, originally denoting the spatial arrangement of cross-sectional units in the sample (Anselin & Bera, 1998). However, the true W is unknown, which is also regarded as the first identification problem in McMillen (2012) and discussed further in the paper of Halleck Vega & Elhorst (2015). Specifically, they criticized that if the pre-specified binary contiguity matrix is taken for granted rather than tested, such pre-specification will lose much infor-mation. There are different ways to specify W, but one methodological step forward would be parameterizing the elements of the W matrix rather than to specify them in advance. Halleck Vega & Elhorst (2015) and Elhorst & Halleck Vega (2017) show that the W matrix within the spatial lag of X (SLX) model, a model accommodating spatial spillover effects generated by the explanatory variables, can easily be parameterized and used to estimate the strength of distance decay effects empirically. At the same time, they mathematically show that this parameterization cannot be used in other spatial econometric models because a perfect solution problem would then occur. This means a global maximum that generates a 100% fit, but which is not of interest from an econometric-theoretical viewpoint (Lee, 2004). Further discussion can be found in Kelejian & Piras (2017).

A lower triangular form of W can be used to solve the perfect solution problem. In the real estate literature, based on the theory that there exists a unidirectional effect of previous housing prices on future prices, Dube & Legros (2016) construct a triangular form of the W matrix that contains two layers of information, spatial and temporal. But a third dimension, similarity, should also be considered. Des Rosiers et al. (2011) and Thanos et al. (2016) stress the importance of ‘comparable sales’, but they only consider the effects of previous similar transactions rather than qualifying and quantifying the similarity. Inspired by the work of Op’t Veld et al. (2008), a measurement of similarity based on degree of convertibility of three groups of housing characteristics, Duran & Elhorst (2018) construct a lower triangular form of W matrix containing three layers of information: spatial, temporal, and similarity. In this paper, we will construct a semi-parameterized three-dimensional W matrix: the strength of distance and time decay effects are parameterized and estimated using non-linear estimation techniques within a SAR model framework, while the degree of similarity with other units is derived from pre-specified information.

(3)

even after the time dimension is added, i.e. the Moran’s I index remains significant after the time adjustment is added. We will use the Moran’s I index to evaluate our parameterized three-dimensional W matrix.

Our paper contributes to literature in several ways. We propose a better W matrix in that it can be parameterized and contains more information. We start with the parameterized matrix 1

dγ_tδ_s, where d represents the distance between the price of the focal house to be

explained and each house that has been sold before in the neighborhood; t the time that has passed between these two transactions; and s the degree of similarity between the two houses. These three layers of information are important since both sellers and buyers will use similar houses that have been sold before within the neighborhood of the focal house as a reference point before making any new bargain. Previous studies only consider transactions occurring within certain distances and time periods, which are usually cut-off values, or they only parameterize distances and time periods1 _{without considering the similarity dimension.}

In this paper, we estimate the parameters γ and δ simultaneously without losing the sim-ilarity information. Our semi-parameterized spatio-temporal-simsim-ilarity W matrix will have a lower triangular form due to the unidirectional effects of the transactions, and this lower triangular form will solve the perfect solution problem. By calculating the Moran’s I index using our semi-parameterized spatio-temporal, we can evaluate the performance of this new W matrix for any remaining spatial dependence.

The remainder of this paper is organized as follows. In Section 2, we review existing liter-ature for different forms of W matrices and Moran’s I index, respectively. In Section 3, we present how to construct our semi-parameterized spatio-temporal-similarity W matrix, how to estimate these two decay parameters, and how to compute the Moran’s I index. Section 4 presents the empirical data, empirical model, and the estimation results. Finally, we provide some further discussions and draw conclusions.

(4)

2 Literature review

Originally, the W matrix is a matrix that reflects the multidirectional spatial relations among the cross-sectional units. Getis & Aldstadt (2004) introduce more than a dozen types of W matrices. The spatial relation is mainly based on nearness principles and distance notions. The simplest W matrix is the binary contiguity matrix (1 for the contiguous neighbors and 0 otherwise) 2 according to the nearness principle. Other matrices borrowing the concept of

nearness are p-order binary contiguity matrices and q-nearest neighbor matrices. W matrices based on the distance notion is inspired by the theory that nearer units have larger effects than those far away, which is known as Tobler’s first law (Tobler, 1970). The most com-mon one is the inverse distance matrix. Other distance decay or decline functions are also proposed, such as bandwidth distance decay (Fotheringham et al., 1996), Gaussian distance decline (LeSage, 2004), and ‘tri-cube’ distance decline (McMillen & McDonald, 2004). The distance can also be parameterized in different ways, such as inverse distance to a universal distance decay power or negative exponential function. Different distance decay parameters for different spatial lags are also allowed (Elhorst & Halleck Vega, 2017). Another spa-tial association is based on descriptive statistics of the data allowing the data to speak for themselves. One example is the spatial cluster using spatial cluster identification techniques (Aldstadt & Getis, 2006). The construction of W matrices could also be driven by economic theories. The most famous one is the gravity function and an application includes the pop-ulation and gross product of each unit (Elhorst & Halleck Vega, 2017). All these matrices reveal the multidirectional cross-sectional dependence which occurs at a specific point in time or is time-invariant.

Dube & Legros (2013) propose a spatio-temporal W matrix that reflects both multidirec-tional spatial dependence and unidirecmultidirec-tional temporal relation given that real estate data is collected over time in the context of large amounts of cross-sectional units. Their empirical results are used to show that the missing of temporal dimension will bias the estimation results (e.g., Dube & Legros, 2014; Dube et al., 2014; Dube et al., 2018). Dube & Legros (2016) justify the lower-triangular form of spatio-temporal W matrix since only past transac-tions have effects on future transactransac-tions while not vice versa in the housing price bargaining

(5)

process. Another property of real estate data is that when modeling housing price inter-actions, only comparable sales should be considered. Based on this argument, Duran & Elhorst (2018) propose a lower-triangular spatio-temporal-similarity W matrix in which a measurement of similarity developed by Op’t Veld et al. (2008) is involved.

The Moran’s I index is one of the most powerful tools for detecting spatial autocorrelation. The global Moran’s I index is applied to test for the spatial autocorrelations, while the local Moran’s I index is calculated to search for spatial clusters (Fu et al., 2011). The calculation of the global Moran’s I index is dependent on the exogenous specification of the W matrix. For a given W matrix, the global Moran’s I index is more stable than other statistics (Florax & Rey, 1995). Dube & Legros (2013) show that the value of the global Moran’s I index becomes smaller suggesting lower spatial autocorrelation when the spatio-temporal W matrix is applied. However, their statistic remains significant suggesting that the spatial autocorrelation is not fully reflected by the spatio-temporal W matrix.

3 Methodology

3.1 Semi-parameterized spatio-temporal-similarity W matrix

To construct our semi-parameterized spatio-temporal-similarity W matrix, we follow the work of Duran & Elhorst (2018). A general N × N parameterized spatio-temporal-similarity W matrix can be split into three parameterized matrices of the same dimensions, where N is the number of the cross-sectional units. The first matrix, W∗_{, denotes the unidirectional}

spatial correlations of the units. The second matrix, T , expresses the unidirectional tem-poral relations among the units. The third matrix, S contains the unidirectional similarity information of the units.

The first step is to construct the parameterized spatial weight matrix, W∗_{. The elements}

on the diagonal are zero since a unit itself will not be considered to be a neighbor. The off-diagonal elements have an inverse distance form to the power of γ larger than 0, which is called the distance decay parameter. Specifically,

(6)

where dij is the distance between unit i and j, γ is the universal distance decay parameter

to be estimated, and ¯d is the distance cut-off value set at 10 kilometres, which is consistent

with the study of Duran & Elhorst (2018). The general form of the parameterized spatial weight matrix W∗ _{is then given by}

W∗=          0 0 0 . . . 0 w∗₂₁ 0 0 . . . 0 ... ... ... ... ... w_{N 1}∗ w_{N 2}∗ w_{N 3}∗ . . . 0          . (2)

It is worth noticing that if γ = 0 the spatial weight matrix is equal to the binary contiguity matrix, and if γ = 1, the spatial weight matrix is defined by a non-parametric inverse dis-tance function.

The next step is to construct the parameterized temporal matrix, T . The diagonal elements are zero duo to the same argument. Dube & Legros (2013) argue that the upper triangular of

T represents the extent to which the past is affected by the future, and the lower triangular

represents the extent to which the past can affect the future when the observations are scaled by time in ascending order. Due to the property of the real estate data, the temporal matrix should have a lower-triangular form. The general elements of the temporal matrix are defined as t∗_ij =      1 tδ ij tij < ¯t 0 tij ≥ ¯t (3) where tij is the time span between unit i and j, δ is the positive time decay parameter to be

estimated, and ¯t is the time period cut-off value set at 13 months, which is inspired by the work of Dube & Legros (2013). They argue that both previous and future transactions have effects on the present selling price. However, when making a purchase decision, buyers tend to compare the selling prices with those of previous transactions. Thus, we only consider the transactions that occurred within the past 13 months. The general form of the temporal matrix T in a lower-triangular form is given by

(7)

The construction of the non-parameterized similarity matrix S needs more work. The off-diagonal elements of S are defined as

s∗_ij =    1 sij sij <¯s 0 sij ≥ ¯s (5) where sij is the similarity index measuring how similar unit i is to unit j. The construction of

the similarity index is inspired by the work of Op’t Veld et al. (2008). The housing character-istics are divided into two groups: primary group and secondary group. Construction area, number of floor, and whether the house is newly built, etc. are part of the primary group. Characteristics such as housing type, whether the house is leased, and has a garage are included in the secondary group. The members of the secondary group are aggregated into one measurement by calculating the mean of each characteristic and computing the pairwise Euclidean distance between the means and each characteristic. Subsequently the aggregated measurement is added to the primary group. The similarity index is then derived by cal-culating the Mahalanobis distance between unit i and j using the updated primary group characteristics. The concept of Mahalanobis distance is brought by Mahalanobis (1936) measuring the distance of a point and the mean of a normal distribution. The advantage of Mahalanobis distance is that it accounts for the correlation in the data by incorporating the variance-covariance matrix of the data (De Maesschalck et al., 2000). Since the similarity is measured by Mahalanobis distance, more similar houses will have smaller Mahalanobis distances to the focal house and therefore will have smaller similarity index values. More similar houses are given more weight because the inverse similarity structure, 1

sij, is larger for

similar houses. The threshold ¯s indicates that only the 50 most similar houses are considered in our case. The final form of the similarity matrix is then given by

S =          0 0 0 . . . 0 s∗₂₁ 0 0 . . . 0 ... ... ... ... ... s∗_{N 1} s∗_{N 2} s∗_{N 3} . . . 0          . (6)

(8)

non-parameterized similarity matrix, which yields a general form, W = W∗ T S          0 0 0 . . . 0 w₂₁∗ t∗₂₁s∗₂₁ 0 0 . . . 0 ... ... ... ... ... w∗_{N 1}t∗_{N 1}s∗_{N 1} w_{N 2}∗ t_{N 2}∗ s∗_{N 2} w∗_{N 3}t∗_{N 3}s∗_{N 3} . . . 0          . (7)

The W matrix with dimension N × N will have a lower-triangular form due to the unidirec-tional relation expressed in the temporal matrix T . Row-normalization is a standard proce-dure in spatial econometrics and matrix normalization is an alternative to row-normalization. The advantage of matrix normalization is that the W matrix will not lose the mutual prop-erties between the elements of the matrix (Elhorst, 2014). We do not apply matrix nor-malization nor row-nornor-malization in this paper for three specific reasons. First, matrix normalization requires dividing each elements of W by its largest eigenvalue. However, for triangular matrices, the eigenvalues are equal to the diagonal elements, which are zero. Sec-ond, traditional row-normalization, in which each element is divided by its row sum, might be an option, but the decay parameters will lose their original economic interpretations. Lastly, the calculation of the Moran’s I index is based on a non-normalized W matrix. If comparable spillover effects measuring the effects of other units on the focal unit are needed, an alternative row-normalization could be used. In particular, the average row sum is set at one. Specifically, the row sums are added together, and the element of the W matrix is rescaled by multiplying the ratio of the sample size to the summation of row sums.

3.2 SAR model

The economic theory behind the cross-sectional spatial autoregressive model (SAR) is ex-plained by the time-dependence motivation of economic agents (LeSage & Pace, 2009, p25). In the long-run equilibrium the SAR model is given by

Y = ρW Y + Xβ + ε, (8)

(9)

the autoregressive parameter (Drukker et al., 2013). The SAR model is usually estimated via maximum likelihood estimation (ML) rather than ordinary least squares (OLS) because OLS estimators are inconsistent due to the endogeneity of W Y , the fact that units affect each other mutually. Halleck Vega & Elhorst (2015) and Elhorst & Halleck Vega (2017) show that the estimation of the SAR model with a parameterized W is far more complicated due to the perfect solution problem, which represents a global maximum that generates a 100% fit, but which is not of interest from an economic viewpoint.

The perfect solution problem can be avoided when our spatial-temporal-similarity W matrix is used due to its lower-triangular form. The SAR model can be used to estimate the two decay parameters via OLS rather than ML becuase the Jacobian term disappears based on the lower-triangular form of the W matrix. Further more, we apply nonlinear estimation techniques to estimate the distance and time decay parameters. Our W matrix also contains the temporal dimension, and thus does not lose since our data set is not purely cross-sectional but has pseudo panel structure. More details will be discussed in the next section. The difference of ML and OLS estimators comes from the Jacobian term in the log-likelihood function of the SAR model, log|I−ρW | (I is the identity matrix). The Jacobian term is equal to zero because of the lower-triangular form of our W matrix. Therefore, ML and OLS will yield the same estimation results. The Hessian matrix used to make statistical inference is calculated by numerical methods instead of analytical methods since the numerical methods take less time and perform better (LeSage & Pace, 2009, p56). A nonlinear estimation technique minimising the residual sum of squares applied by Halleck Vega & Elhorst (2015) is modified to estimate those two decay parameters.

3.3 Global Moran’s I index

(10)

I index. Following Dube & Legros (2013), the global Moran’s I index and relevant statistics can be derived as follows,

M(I) = N W0

r0W r

r0r (9)

where W0 is the sum of elements of the non-normalized W matrix and r 0

is the residual vector. The mean (E(I)) and variance (V ar(I)) of this index are given as follows,

E(I) = tr(MW ) N − K V ar(I) = tr(MW MW 0 ) + tr(MW )2_{+ [tr(MW )]}2 (N − K)(N − K + 2) , (10) where M is equal to I−PX and PX is a projection matrix. M(I) has a normal distribution

with mean E(I) and variance V ar(I). Statistical inference can be drawn from the statistic given by

z = Mq(I) − E(I)

V ar(I)

. (11)

Under the null hypothesis, the global Moran’s I index is invariant when the hedonic model is estimated by OLS (King, 1981). The null hypothesis states that there exists no spatial correlation. Within the SAR framework, this statistic should be insignificant after controlling for spatial dependence.

4 Empirical analysis

4.1 Data

The housing transaction data include selling price per square metre, sales date, geographical coordinates, and housing characteristics spanning from 1993 to 2017 in the Dutch city of Groningen and are from the Dutch Association of Realtors (NVM)3. Some houses might be

sold several times, but repeated sales only account for a very small part of all the transac-tions and can be ignored. Therefore, the structure of our data set takes the form of a pseudo panel. The data of one particular type of house, half-of-double, which accounts for 5.64% of all sales are applied to estimate our semi-parameterized three-dimensional W matrix. In total, our data set covers 2722 transactions conducted in 29 neighbourhoods (distinguished

(11)

by postal codes) in the city of Groningen from 1993 to 2017. Data sets of two other types of house, corner and detached, cannot be used to estimate the distance and time decay pa-rameters. More details will be obtained in the simulation part.

The housing characteristics covered in our dataset and applied in our approach are reported in Table 1. Among them, there are 36 variables, which include two discrete variables, 30 dummy variables, and four continuous variables. Discrete variables are number of rooms and floors. Continuous variables are construction area, land area, garden area, and space size. Discrete variables include attic presence, location types, construction periods, etc.

Description Type Description Type

Construction area (in m2) Continuous Presence of balcony Dummy Land area (in m2₎ _Continuous _{Whether the house is leased} _Dummy

Garden area (in m2) Continuous Presence of attic Dummy Space size (in m3) Continuous Presence of cellar Dummy Number of floors Discrete Whether it is a monument Dummy Number of rooms Discrete Whether it is newly built Dummy

Location type 1 Dummy Maintenance(inside 5-6) Dummy

Location type 2 Dummy Maintenance(inside 7+) Dummy

Location type 3 Dummy Maintenance(outside 5-6) Dummy

Location type 4 Dummy Maintenance(outside 7+) Dummy

Construction period (1906-1930) Dummy Presence of isolation Dummy Construction period (1931-1944) Dummy Presence of gas or coal stove Dummy Construction period (1945-1959) Dummy Presence of boiler Dummy Construction period (1960-1970) Dummy Presence of solar Dummy Construction period (1971-1980) Dummy Presence of fireplace Dummy Construction period (1981-1990) Dummy Presence of garage Dummy Construction period (1991-2000) Dummy Presence of barn Dummy Construction period (after 2001) Dummy Presence of parking facility Dummy

Table 1: Descriptions of Housing Characteristics

(12)

This trend suggests that year dummies should enter the model to control for business cycle effects. Later, the statistical inference results also confirm the validity of year dummies.

Figure 1: Annual Average Selling Price Per Square Metre

4.2 Empirical models

The empirical model within the SAR model framework is given by

log(SPit) = ρ n

X

j=1

w_ij∗t∗_ijs∗_ijlog(SPjt) + Xitβ+ µr+ λt+ εit. (12)

The dependent variable is the selling price per square metre in log form for observation i at time t. Xit is a row vector containing the control variables 4 reported in Table 1 for

observation i at time t. µr and λt are neighbourhood and year dummies to control the

neighbourhood and time fixed effects, respectively. εit is the iid error. N is the number of

the observations.

(13)

4.3 Simulation results

To roughly test whether we can obtain a maxima of the log-likelihood value, we conduct simulation studies by calculating the log-likelihood values when setting different values of distance (time) decay parameter while keeping non-parameterized time (distance) and simi-larity matrices unchanged. Both of the parameters vary from 0 to 10 with increments equal to 0.1. Figures 2 to 4 depict the log-likelihood values with respect to different pre-set values of distance and time decay parameters for different types of houses. For the corner house type, the log-likelihood values with respect to the distance decay parameter have a maxima with a Kinky point, while the log-likelihood values with respect to the time decay parameter do not have a clear pattern. For the corner house type, neither of the curves is concave. As for the half-of-double house type, both of the curves are concave after certain points even though the points of the log-likelihood values with respect to distance decay parameter are rather flat after reaching the maxima. With the simulation results, we can infer that only the data of the half-of-double houses can be applied to parameterization. One explaination could be that only the data of the city of Groningen are used and the housing prices of other regions could also have an effect. Therefore, the data of a larger area might be used. Another explaination could be that the specifications of the parameterized spatial weight matrix and temporal matrix are baesd on inverse ‘distance’ function. Other specifications, such as negative exponential fuction, can be an alternative. Then we proceed to formally estimate the distance and time decay parameters.

(14)

Figure 3: Simulation Results of House Type Detached

Figure 4: Simulation Results of House Type Half of Double

4.4 Estimation Results

Table 2 shows the estimation of the autoregressive parameters, the distance, and time decay parameters and corresponding t-statistics 5_{. Column (1) reports the results when neither}

neighbourhood (nbh) nor time-period effects are added; columns (2) and (3) report the re-sults when either neighbourhood or time-period fixed effects are included in the model, and

(15)

column (4) 6 _{reports the results when both neighbourhood and time-period effects are}

in-cluded. For models without fixed effects or only neighbourhood effects, only the time decay effect is significant at the 5% significance level. The estimated decay parameters become much larger for the models with only time-period and both fixed effects, and both of the pa-rameters are significant at the 5% significance level for these two models. All the estimated autoregressive parameters are significant at 5% significance level. The value of the estimated autoregressive parameter is positive and small when both effects are included suggesting that similar houses will generate positive spatial-temporal spillover effects. The LR tests, which are equal to two times the value of the difference between the log-likelihood value of the un-restricted model and that of the un-restricted model, show that both effects should be included. The statistics of the LR test have a chi-squared distribution, and the degrees of freedom are equal to the number of the restrictions that need to be imposed. The null hypothesis that the neighbourhood or year dummies are jointly insignificant can be rejected, given that the LR test results are 972.9 with 25 degrees of freedom and 1025.6 with 29 degrees of freedom, respectively. Consequently, our estimated distance and time decay parameters are 3.821 and 3.294 while keeping the pre-specified similarity matrix unchanged. Similar houses located far away or evaluated long ago receive lower weights. R2 _{statistics are not reported here}

since R2 _{cannot be used to select the W matrix when the specification of the W matrix is}

purely empirical analysis with uncertainty (Kelejian & Piras, 2017, p90).

(1) (2) (3) (4)

Est. T-stat. Est. T-stat. Est. T-stat. Est. T-stat.

ρ -0.320 -2.831 -0.290 -3.606 0.000 7.651 0.00006 5.493

γ 0.000 0.000 0.080 0.706 4.485 71.649 3.821 25.565

δ 1.642 11.946 1.616 16.085 4.010 17.245 3.294 7.191

Nbh fixed effects No Yes No Yes

Time-period fixed effects No No Yes Yes

Log-likelihood -3121.3 -2125.1 -2177.8 -1152.2

Table 2: Estimation Results for Distance and Time Decay Parameters

(16)

After constructing our semi-parameterized three-dimensional W matrix, we need to evaluate our W matrix. The techniques proposed by Dube & Legros (2013) are applied in our paper. The residuals are calculated based on the hedonic price model, which is the regression of housing price (per square metre) in log form on spatial autocorrelation terms and housing characteristics listed above. Moran’s I indices calculated using different forms of W matri-ces are reported in Table 3. These statistics have a normal distribution when the hedonic price model is estimated by OLS. All three indices are insignificant at the 5% significance level, suggesting that there is no spatial correlation using these three W matrices. The non-parameterized three-dimensional W matrix has a lower Moran’s I index value than that of the spatio-temporal W matrix, suggesting that without the similarity information, the specifica-tion of the model might be inaccurate. The Moran’s I index based on our semi-parameterized three-dimensional W matrix has the lowest value, declaring that the estimation model is the most accurate if the newly proposed W matrix is applied. Besides, our proposed W matrix is one step further than the pure spatio-temporal and non-parameterized three-dimensional W matrices since it suffers less from the pre-specification problem.

Moran’s I index Z-score

Spatio-temporal W matrix 0.1409 0.00018

Spatio-temporal-similarity W matrix 0.1213 0.0014

Semi-parameterized spatio-temporal-similarity W matrix 0.0260 0.00007 Table 3: Moran’s I index

5 Conclusion and Discussion

(17)

econometric studies (Elhorst, 2014) since the interpretation of spillover effects are based on it. Row-normalization with the average row sum set equal to one can work as an alternative to the matrix normalization and the traditional row-normalization.

(18)

References

Aldstadt, J., & Getis, A. (2006). Using amoeba to create a spatial weights matrix and identify spatial clusters. Geographical analysis, 38 (4), 327–343.

Anselin, L., & Bera, A. K. (1998). Introduction to spatial econometrics. Handbook of Applied Economic Statistics,237.

De Maesschalck, R., Jouan-Rimbaud, D., & Massart, D. (2000). The mahalanobis distance.

Chemometrics and intelligent laboratory systems, 50 (1), 1–18.

Des Rosiers, F., Dube, J., & Theriault, M. (2011). Do peer effects shape property values?

Journal of Property Investment & Finance, 29 (4/5), 510–528.

Drukker, D., Egger, P., & Prucha, I. (2013). On two-step estimation of a spatial autoregres-sive model with autoregresautoregres-sive disturbances and endogenous regressors. 32 (5-6).

Dube, J., & Legros, D. (2013). A spatio-temporal measure of spatial dependence: An example using real estate data. Papers in Regional Science, 92 (1), 19–30.

Dube, J., & Legros, D. (2014). Spatial econometrics and the hedonic pricing model: what about the temporal dimension? Journal of Property Research, 31 (4), 333–359.

Dube, J., & Legros, D. (2016). A spatiotemporal solution for the simultaneous sale price and time-on-the-market problem. Real Estate Economics, 44 (4), 846–877.

Dube, J., Legros, D., & Thanos, S. (2018). Past price ‘memory’ in the housing mar-ket: testing the performance of different spatio-temporal specifications. Spatial Economic

Analysis, 13 (1), 118–138.

Dube, J., Legros, D., Theriault, M., & Des Rosiers, F. (2014). A spatial difference-in-differences estimator to evaluate the effect of change in public mass transit systems on house prices. Transportation Research Part B: Methodological, 64 , 24–40.

Duran, N., & Elhorst, J. P. (2018). A spatio-temporal-similarity and common factor approach of individual housing prices.

(19)

Elhorst, J., & Halleck Vega, S. (2017). The slx model: Extensions and the sensitivity of spatial spillovers to w.

Florax, R., & Rey, S. (1995). The impacts of misspecified spatial interaction in linear

regres-sion models. Springer, Berlin, Heidelberg.

Fotheringham, A., Charlton, M., & Brunsdon, C. (1996). The geography of parameter space: an investigation of spatial non-stationarity. International Journal of Geographical

Information Systems, 10 (5), 605–627.

Fu, W., Zhao, K., Zhang, C., & Tunney, H. (2011). Using moran’s i and geostatistics to identify spatial patterns of soil nutrients in two different long-term phosphorus-application plots. Journal of Plant Nutrition and Soil Science, 174 (5), 785–798.

Getis, A., & Aldstadt, J. (2004). Constructing the spatial weights matrix using a local statistic. Geographical analysis, 36 (2), 90–104.

Halleck Vega, S., & Elhorst, J. (2015). The slx model. Journal of Regional Science, 55 (3), 339–363.

Kelejian, H., & Piras, G. (2017). Spatial econometrics. Academic Press.

King, M. (1981). A small sample property of the cliff-ord test for spatial correlation. Journal

of the Royal Statistical Society: Series B (Methodological), 43 (2).

Lee, L. (2004). Asymptotic distributions of quasi-maximum likelihood estimators for spatial autoregressive models. Econometrica, 72 (6), 1899–1925.

LeSage, J. (2004). A family of geographically weighted regression models. Springer. LeSage, J., & Pace, R. (2009). Introduction to spatial econometrics. CRC Press.

Mahalanobis, P. (1936). , 1936. on the generalized distance in statistics. National Institute

of Science of India.

McMillen, D., & McDonald, J. (2004). Locally weighted maximum likelihood estimation:

Monte Carlo evidence and an application. Springer.

(20)

Op’t Veld, D., Bijlsma, E., & van de Hoef, P. (2008). Automated valuation in the dutch housing market: The web-application ‘marktpositie’ used by nvm-realtors. Advances in

Mass Appraisal Methods, 58 , 70–90.

Thanos, S., Dube, J., & Legros, D. (2016). Putting time into space: the temporal coherence of spatial applications in the housing market. Regional Science and Urban Economics, 58 , 78–88.

Tobler, W. R. (1970). A computer movie simulating urban growth in the detroit region.

(21)

A

Appendix: Estimation Results with Two Fixed

Ef-fects

Est. T-stat. ρ 0.00006 5.493 M2 -0.004 -12.819 n-stores 0.079 7.060 lot area 0.0002 6.347 n-rooms 0.003 0.581 n-balc 0.071 4.706 loc garden 0.000 -0.538 lease -0.006 -0.064 M3 0.0006 8.434 ATTIC -0.025 -2.357 cellar 0.134 4.187 MONUMENT 0.562 9.535 newly built -0.155 -3.153 GARAGE 0.033 2.142 barn 0.006 0.553 no parking facility 0.011 0.609 location beautiful 0.021 1.963 location busy road -0.053 -1.786

location in centre 0.506 11.510 location in residential area -0.015 -1.117

(22)

>2001 0.273 6.975 Maintenance inside 5-6 0.523 7.272 Maintenance inside 7+ 0.641 8.874 Maintenance outside 5-6 0.869 11.630 Maintenance outside 7+ 0.887 11.782 ISOLATION -0.012 -0.741

Gas or coal stove 0.496 9.966

(23)

(24)

B

Appendix: Matlab Codes

% CONSTRUCTION OF THREE MATRICES

A = xlsread ('X:\ My Desktop \ half \ half . xlsx ');

% vif(A) % Variance inflation factor

nobs = length (A);

Check = zeros (nobs ,1);

% M2 values in the brackets are edge values

A1= discretize (A(: ,4) ,[0 50 100 150 200 1000]) ;

A2= discretize (A(: ,15) ,[0 50 100 200 300 400 500 1000 100000]) ;

% M2 land

A3= discretize (A(: ,18) ,[0 50 100 200 300 400 500 1000 1000000]) ;

% M2 garden

A4= discretize (A(: ,21) ,[0 100 200 300 400 500 600 1000 1500 10000]) ; % M3

A(: ,4)=A1; A(: ,15)=A2; A(: ,18)=A3; A(: ,21)=A4;

% Note original A is out of memory now !!!! % column 20 is deleted due to all zeros

col1 =[4 14 15 16 18 21 22 23 26];

% column 24 is deleted due to all zeros % column 12 is deleted

col2 =[19 25 27 28 29 30 31 32 33 47 48 49 50 51 52 53 54 55]; MEAN = mean (A(:, col2 ));

X2dist = pdist2 (MEAN ,A(:, col2 )) '; Xsim =[A(:, col1 ) X2dist ];

COV=cov( Xsim );

(25)

% Construction of W matrices % % sorttype =A(: ,12); sortdate =A(: ,8); sortx1 =A(: ,9); sortx2 =A(: ,10);

Wdist = speye ( nobs ); % identity matrix

Wtime = speye ( nobs ); Wsim = speye ( nobs ); i=1;

% maximum time length is set at one year to limit W

while sortdate (i) >403 %&& i <10

count =0; j=i+1;

while sortdate (i)-sortdate (j)==0

j=j+1; end

Jstart =j;

while sortdate (i)-sortdate (j) <403 && j< nobs % % Dube study

1000 m and 13 months

distance = sqrt (( sortx1 (i)-sortx1 (j)) ˆ2+( sortx2 (i)-sortx2 (j))ˆ2) /1000;

if ( distance <10) % houses should be of the same type .

count = count +1;

Wdist (i,j)=1/ distance ;

Wtime (i,j) =1/( sortdate (i)-sortdate (j));

Wsim (i,j)=1/ pdist2 ( Xsim (j ,:) ,Xsim (i ,:) ,'mahalanobis ',COV);

if Wdist (i,j) >1000 Wdist (i,j) =1000; end % some values are infinities

(26)

j=j+1; end Jend =j; Check (i ,1)= count ; i=i+1; end

Wdist =Wdist - speye ( nobs ); % Matrix is not normalized yet

% Wdist = Wdist ./ sum(Wdist ,2); % Now it is normalized , but this might not be necessary

Wtime =Wtime - speye ( nobs ); Wsim =Wsim - speye ( nobs );

save ('X:\ My Desktop \ half \ Wdist .mat ','Wdist ') save ('X:\ My Desktop \ half \ Wtime .mat ','Wtime ') save ('X:\ My Desktop \ half \ Wsim .mat ','Wsim ')

% SELECT THE 50 MOST SIMILAR HOUSES

load ('X:\ My Desktop \ half \ half \ Wsim ') Wsim_50 = Wsim ;

N = length ( Wsim ); n2 = N -50;

[b ind ]= sort (Wsim ,2); for i = 1:N

for j = 1:N

if ind(i,j)>n2

Wsim_50 (i,j) = Wsim (i,j); else

Wsim_50 (i,j)= 0; end

(27)

save ('X:\ My Desktop \ half \ Wsim_50 .mat ','Wsim_50 ')

% SIMULATION RESULTS

B = xlsread ('X:\ My Desktop \ half \ half . xlsx '); load ('Wdist ')

load ('Wtime ') load ('Wsim_50 ')

% generate regional and time dummies

nobs = length (B); N = nobs ;

(28)

year15 = zeros (nobs ,1); year16 = zeros (nobs ,1); year17 = zeros (nobs ,1); region12 = zeros (nobs ,1); region13 = zeros (nobs ,1); region14 = zeros (nobs ,1); region15 = zeros (nobs ,1); region16 = zeros (nobs ,1); region17 = zeros (nobs ,1); region18 = zeros (nobs ,1); region21 = zeros (nobs ,1); region22 = zeros (nobs ,1); region23 = zeros (nobs ,1); region24 = zeros (nobs ,1); region25 = zeros (nobs ,1); region26 = zeros (nobs ,1); region27 = zeros (nobs ,1); region28 = zeros (nobs ,1); region31 = zeros (nobs ,1); region32 = zeros (nobs ,1);

% region33 = zeros (nobs ,1);

(29)

for i=1: nobs

(30)

if (B(i ,2) ==9722) region22 (i) = 1; end if (B(i ,2) ==9723) region23 (i) = 1; end if (B(i ,2) ==9724) region24 (i) = 1; end if (B(i ,2) ==9725) region25 (i) = 1; end if (B(i ,2) ==9726) region26 (i) = 1; end if (B(i ,2) ==9727) region27 (i) = 1; end if (B(i ,2) ==9728) region28 (i) = 1; end if (B(i ,2) ==9731) region31 (i) = 1; end if (B(i ,2) ==9732) region32 (i) = 1; end

%if (B(i ,2) ==9733) region33 (i) = 1; end

if (B(i ,2) ==9734) region34 (i) = 1; end if (B(i ,2) ==9735) region35 (i) = 1; end if (B(i ,2) ==9736) region36 (i) = 1; end if (B(i ,2) ==9737) region37 (i) = 1; end if (B(i ,2) ==9738) region38 (i) = 1; end if (B(i ,2) ==9741) region41 (i) = 1; end if (B(i ,2) ==9742) region42 (i) = 1; end if (B(i ,2) ==9743) region43 (i) = 1; end if (B(i ,2) ==9744) region44 (i) = 1; end if (B(i ,2) ==9745) region45 (i) = 1; end if (B(i ,2) ==9746) region46 (i) = 1; end end

region_fix = [ region12 region13 region14 region15 region16 region17 region18 ...

region21 region22 region23 region24 region25 region26 region27 ...

region28 region31 region32 region34 region35 region36 region37 region38 ...

region41 region42 region43 region44 region45 region46 ];

(31)

year03 year04 year05 year06 year07 year08 year09 year10 year11 year12...

year13 year14 year15 year16 year17 ]; [ rfe_row rfe_col ]= size ( region_fix ); [ tfe_row tfe_col ]= size ( time_fix ); B_r = [B region_fix ];

B_t = [B time_fix ];

B_both = [B region_fix time_fix ];

Y = log(B(1: end ,5));% average price per m2

x1=B(1: end ,4); [b1 b2 ]= size (B);

[ b_r_1 b_r_2 ]= size (B_r); [ b_t_1 b_t_2 ]= size (B_t);

[ b_both_1 b_both_2 ]= size ( B_both ); x2= B_both (1: end ,14:19) ;

x3= B_both (1: end ,21:23) ; x4= B_both (1: end ,25:41) ;

x5= B_both (1: end ,47: b_both_2 );

X_both = [x1 x2 x3 x4 x5 ];% should include a column vector of ones , plz check the ols fuction

[ both_row both_col ]= size ( X_both ); X = X_both ;

count =0;

%for p =5.5

for p =0.0:0.1:10 count = count +1;

W2 =( Wdist .ˆp).* Wtime .* Wsim_50 ;

%W2=W2/max( eigs (W2));

(32)

results =ols(Y ,[ WY X]);

lik =-(N/2)*log (2* pi) -(N/2)*log( results .si2) -(1/ results .si2)* results .resid '* results . resid ;

score (count ,:) =[p lik ]; end

x = score (1: end ,1); y = score (1: end ,2); fig = plot (x,y);

saveas (fig , 'distance_50 .pdf ')

clear count score y x p lik resid si2 WY fig count =0;

%for p =5.5

for p =0.0:0.1:10 count = count +1;

W2 =( Wtime .ˆp).* Wdist .* Wsim_50 ;

%W2=W2/max( eigs (W2));

WY = W2*Y;

results =ols(Y ,[ WY X]);

lik =-(N/2)*log (2* pi) -(N/2)*log( results .si2) -(1/ results .si2)* results .resid '* results . resid ;

score (count ,:) =[p lik ]; end

x = score (1: end ,1); y = score (1: end ,2); fig = plot (x,y);

(33)

clear count score y x p lik resid si2 WY fig % ESTIMATION RESULTS x1=B(1: end ,4); [b1 b2 ]= size (B); [ b_r_1 b_r_2 ]= size (B_r); [ b_t_1 b_t_2 ]= size (B_t);

[ b_both_1 b_both_2 ]= size ( B_both ); x2=B(1: end ,14:19) ;

x3=B(1: end ,21:23) ; x4=B(1: end ,25:41) ; format long g

% %%%%%%%%%%%%%%%%%%%%%%%%%% % neither region nor time % % %%%%%%%%%%%%%%%%%%%%%%%%%%

x5=B(1: end ,47: b2);

X_non = [x1 x2 x3 x4 x5 ];% should include a column vector of ones , plz check the ols fuction

[ non_row non_col ]= size ( X_non ); X = X_non ;

iter =0; converge =1.0; criteria =0.0001; itermax =25; options . Display ='off ';

(34)

[alpha ,fval , exitflag ]= fminsearch ('only_two ',alpha0 , options ,Y,X, Wdist ,Wtime , Wsim );

gam1_non = alpha (1); gam2_non = alpha (2);

W_nonnor_non = ( Wdist .ˆ gam1_non ).*( Wtime .ˆ gam2_non ).* Wsim ; WY = W_nonnor_non *Y;

X_all = [WY X];

results =ols(Y, X_all );

%% prt_reg ( results );

loglikelihood_non =-(N/2)*log (2* pi) -(N/2)*log( results .si2) -(1/ results .si2)* results .resid '* results . resid ;

parm =[ results . beta ;alpha '];

xpx = hessian ('only_two_2 ',parm ,Y,X,Wdist ,Wtime , Wsim ); [ ndummy nvar ]= size ( X_all );

xpx (1: nvar ,1: nvar )=X_all '* X_all ; %(4: x1 , x2 , wx1 , wx2) replace by x'x.

varcov = results . sige * invpd (xpx); % xpx is al positive definite ( sige : e'e/(n-k))

tvalues = parm ./ sqrt (abs( diag ( varcov ))); chart_non =[ parm tvalues ];

clear X X_all parm tvalues alpha WY W_nor WY_nor results ;

% %%%%%%%%%%%%%%%%%%%%%%%%%% % only time fixed effects % % %%%%%%%%%%%%%%%%%%%%%%%%%%

x6=B_t (1: end ,47: b_t_2 );

X_t = [x1 x2 x3 x4 x6 ];% should include a column vector of ones , plz check the ols fuction

(35)

options . MaxFunEvals =1000; options . MaxIter =1000; options . TolX =0.0001; options . TolFun =0.0001;

[alpha ,fval , exitflag ] = fminsearch ('only_two ',alpha0 , options ,Y, X,Wdist ,Wtime , Wsim );

gam1_t = alpha (1); gam2_t = alpha (2);

W_nonnor_t = ( Wdist .ˆ gam1_t ).*( Wtime .ˆ gam2_t ).* Wsim ; WY = W_nonnor_t *Y;

X_all = [WY X];

loglikelihood_t =-(N/2)*log (2* pi) -(N/2)*log( results .si2) -(1/ results .si2)* results .resid '* results . resid ;

tvalues = parm ./ sqrt (abs( diag ( varcov ))); chart_t =[ parm tvalues ];

(36)

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % only regional fixed effects % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x7=B_r (1: end ,47: b_r_2 );

X_r = [x1 x2 x3 x4 x7 ];% should include a column vector of ones , plz check the ols fuction

[ r_row r_col ]= size (X_r); X = X_r;

gam1_r = alpha (1); gam2_r = alpha (2);

W_nonnor_r = ( Wdist .ˆ gam1_r ).*( Wtime .ˆ gam2_r ).* Wsim ; WY = W_nonnor_r *Y;

X_all = [WY X];

loglikelihood_r =-(N/2)*log (2* pi) -(N/2)*log( results .si2) -(1/ results .si2)* results .resid '* results . resid ;

(37)

by x'x.

tvalues = parm ./ sqrt (abs( diag ( varcov ))); chart_r =[ parm tvalues ];

clear X X_all parm tvalues alpha WY W_nor WY_nor results ;

% %%%%%%%%%%%%%%% % both effects % % %%%%%%%%%%%%%%%

gam1_both = alpha (1); gam2_both = alpha (2);

W_nonnor_both = ( Wdist .ˆ gam1_both ).*( Wtime .ˆ gam2_both ).* Wsim ; WY = W_nonnor_both *Y;

X_all = [WY X];

(38)

loglikelihood_both =-(N/2)*log (2* pi) -(N/2)*log( results .si2) -(1/ results .si2)* results .resid '* results . resid ;

tvalues = parm ./ sqrt (abs( diag ( varcov ))); chart_both =[ parm tvalues ];

% %%%%%%%%% % LR test % % %%%%%%%%%

LR_r = 2*( loglikelihood_r - loglikelihood_non ); % test for regional fixed effects

LR_t = 2*( loglikelihood_t - loglikelihood_non ); % test for time fixed effects

LR_both = 2*( loglikelihood_both - loglikelihood_non ); % test for both effects

save ('X:\ My Desktop \ half \ chart_non .mat ','chart_non ') save ('X:\ My Desktop \ half \ chart_r .mat ','chart_r ') save ('X:\ My Desktop \ half \ chart_t .mat ','chart_t ')

save ('X:\ My Desktop \ half \ chart_both .mat ','chart_both ') save ('X:\ My Desktop \ half \ loglikelihood_non .mat ','

loglikelihood_non ')

save ('X:\ My Desktop \ half \ loglikelihood_t .mat ',' loglikelihood_t '

)

(39)

)

save ('X:\ My Desktop \ half \ loglikelihood_both .mat ','

loglikelihood_both ') %MORAN 'S I INDEX x1=B(1: end ,4); [b1 b2 ]= size (B); [ b_r_1 b_r_2 ]= size (B_r); [ b_t_1 b_t_2 ]= size (B_t);

[ b_both_1 b_both_2 ]= size ( B_both ); x2= B_both (1: end ,14:19) ;

x3= B_both (1: end ,21:23) ; x4= B_both (1: end ,25:41) ;

W = ( Wdist .ˆ alpha1 ).*( Wtime .ˆ alpha2 ).* Wsim_50 ; results = ols(Y ,[W*Y X]);

s_0 = sum(sum(W)); %W is not normalized

[N K]= size (X);

e = results . resid ; % residuals of the hedonic price model

M_I = (N/s_0)*(e '*W*e/(e '*e)); M = speye (N) - X*inv(X '*X)*X '; E_M_I = trace (M*W)/(N-K);

Var_M_I =( trace (M*W*M*W ')+ trace ((M*W)ˆ2) +( trace (M*W))ˆ2) /(N-K) *(N-K+2) -E_M_I ˆ2;

(40)

z_1 = z;

save ('X:\ My Desktop \ half \ M_I_1 .mat ','M_I_1 ') save ('X:\ My Desktop \ half \z_1.mat ','z_1 ') clear W M_I M E_M_I Var_M_I z results W = Wdist .* Wtime ;

results = ols(Y ,[W*Y X]);

[N K]= size (X);

M_I = (N/s_0)*(e '*W*e/(e '*e)); M = speye (N) - X*inv(X '*X)*X '; E_M_I = trace (M*W)/(N-K);

z = (M_I - E_M_I )/ Var_M_I ˆ0.5; M_I_2 = M_I;

z_2 = z;

save ('X:\ My Desktop \ half \ M_I_2 .mat ','M_I_2 ') save ('X:\ My Desktop \ half \z_2.mat ','z_2 ') clear W M_I M E_M_I Var_M_I z results W = Wdist .* Wtime .* Wsim_50 ;

results = ols(Y ,[W*Y X]);

[N K]= size (X);

(41)

M = speye (N) - X*inv(X '*X)*X '; E_M_I = trace (M*W)/(N-K);

z = (M_I - E_M_I )/ Var_M_I ˆ0.5; M_I_3 = M_I;

z_3 = z;