The importance of measuring residential property price inflation is paramount to households and economies; however, the heterogeneity of properties makes it difficult (De Haan & Diewert 2011). Log linear models have been used extensively in real estate economics to estimate property prices and inflation. This study investigates generalised linear models as an alternative to the typical log linear approach. Cross-sectional hedonic price functions are developed and compared for the South African residential property market over a 5-year period.

Residential property is an important component of individual and national wealth where it is capitalised on household balance sheets, informing economic policy formulation (Hill 2013). Goodhart and Hoffman (2008) conducted a study providing evidence of relationships between house prices, credit and broad money. Using vector auto-regression fitted with ordinary least squares, their research showed statistically significant relationships between home prices and the macro economy. Bordo and Jeanne (2002) found an increased likelihood of a financial crisis occurring when real estate prices reached a peak or shortly after a bust, in a study of advanced economies spanning from 1970 to 2001. This resonates with the views of De Haan and Diewert (2011) who assert that sharp declines in home prices can adversely affect the debt to equity ratio and credit ratings. Residential property has an important role in economies and understanding price inflation is imperative; however, measuring price inflation is difficult because of infrequent transactions and the heterogeneity of properties.

Hedonic regression is ubiquitous in the construction of residential property price indices where log linear models are commonly developed in the price estimation procedure (De Haan & Diewert 2011; Jiang et al. 2015). Hedonic regression has been found useful as a quality-adjusted methodology where pure price changes are measured and not simply changes in the composition of samples in different periods (Shimizu, Nishimura & Watanabe 2010). Hedonic pricing measures the price of an item through its utility bearing characteristics where the price of the item is determined by the vector of its characteristics (Rosen 1974). Hedonic pricing describes the functional relationship of a heterogeneous item and the implicit attributes: Pj=P(Zj)=P(Zj1,Zj2,…,Zjk)[Eqn 1] where P_j is the price of the jth item which is a function of a set of characteristics Z_j (Goodman 1978). Hedonic pricing is useful when estimating the price of heterogeneous goods. Heterogeneous or differentiated goods are goods that differ in their respective composition of characteristics; however, consumers consider the set of characteristics closely related, defining it as a single item (Day 2003). Hedonic pricing mathematically models residential property prices as a function of structural and location characteristics (Lyons 2015). Ordinary least squares models are typically employed to estimate the marginal contributions of each characteristic, taking the form: θ(μ)=Xβ[Eqn 2] where β is p < n unknown parameters, the matrix X_n×p is a set of known independent variables and Xβ is the linear structure (Lindsey 2005). The implicit price for characteristic i of property j is calculated by taking the partial derivative. Because of the positive domain and positively skewed nature of residential property prices, log linear models are often adopted as a quality-adjusted technique that controls for changes in the quality of properties transacted in different periods, whilst reducing heteroscedasticity in the residuals (Silver 2016). Day (2003); Bourassa, Cantoni and Hoelis (2007); and Els and Von Fintel (2010) conducted separate hedonic price studies for different property markets using log linear models. Els and Von Fintel (2010) found that the assumption of the linear functional form was violated, finding quantile regression more appropriate to capture the hedonic price function. A potential problem with transforming property prices to the log scale is that exponentiation of the fitted values produces geometric mean estimates and not arithmetic mean estimates (Olivier, Johnson & Marshall 2008). Another potential concern is the assumption that property prices are lognormal when a different distribution family may represent the data better. Generalised linear models incorporate exponential classes of distribution families, which facilitate modelling the response on the original scale. The objective of this study is to investigate generalised linear models as an alternative framework to the log linear model in the development of hedonic price functions for the South African residential property market.

Generalised linear models are a natural extension of classical linear models where properties such as linearity and computing parameter estimates are similar (McCullagh & Nelder 1989). Generalised linear models are characterised by three components. Firstly, a stochastic or random component representing a response variable y, consisting of independent observations (y₁, y₂, … , y_n), belonging to a class of an exponential family distribution in the form of: f(y;θ∅)=exp{θy−b(θ)∅+c(y,∅)}[Eqn 3] where Ø is a dispersion parameter and b(.), c(.) are known functions and the range of Y does not depend on θ or Ø. For a random response variable Y with distribution of form 3 E(y) = μ. Secondly, a systematic component that consists of a set of covariates (x₁, x₂, … , x_p) which combine linearly with the coefficients to produce the linear predictor η. Therefore, η = βX. Finally, a link function that connects the stochastic and systematic components where η = µ.

This generalisation takes the form: ηi=g(μi)[Eqn 4] where g (.) denotes the link function and η = μ through the link function. The link function relates the conditional mean to the systematic component, namely the covariates. This formulation allows for the exponential family of distributions including normal; however, the link function may become any monotonic differentiable function, which then allows extensions to distributions such as Poisson, binomial and gamma amongst others (McCullagh & Nelder 1989). This means that generalised linear models are suitable for modelling continuous data as well as count and binary data.

Generalised linear models obtain maximum likelihood estimates of parameters belonging to an exponential distribution family using the iterative reweighted least squared algorithm where the link function makes the systematic effects linear (Nelder & Wedderburn 1972). Maximum likelihood estimates are a vector of parameter estimates produced by a model function which makes the observed data probable given the model function (Lindsey 2005).

The primary goodness-of-fit measure for generalised linear models is called the deviance which is the logarithm of a ratio of likelihoods (McCullagh & Nelder 1989). The analysis of deviance makes model assessment and comparison possible in terms of the choice of covariates. Given a set of data, two extreme models are possible. Firstly, a null model with one parameter which represents a common μ for all the ys. Secondly, a complete model where all the ys are different, matching the data completely. Fitting a model with more than one parameter represents a saturated model that can be compared to the null model (Dobson & Barnett 2018). The fitting of n parameters is performed by maximising the likelihood of matching the model to the likelihood of the data through the deviance that differs based on the distribution.

For the normal distribution, the deviance is simply the sum of squares just like ordinary least squares which means that fitting a normal or lognormal distribution with the identity link function, where the natural logarithm of the response is taken, is equivalent to fitting a linear or log linear ordinary least squares model. For generalised linear models, the saturated model should have a lower deviance than the null model, indicating that the inclusion of n parameters is a better fit. Guisan and Zimmernam (2000) propose that variance reduction in model formulation is generally a desired characteristic of the goodness of fit as with generalised linear models, where deviance reduction can be converted to an equivalent R² statistic: D2=(Null deviance−Residual deviance)/Null deviance[Eqn 5] where D² is the deviance explained or the amount of deviance accounted for by the model. Naturally, this leads to an understanding of the residuals of generalised linear models where the deviance residuals are reported as a measure of discrepancy. Deviance residuals are calculated as follows: sign(yi−μ^)di2[Eqn 6]

This formulation shows that deviance residuals are calculated by taking the signed square root of the ith observation to the total model deviance. One can begin to understand the quality of fit that reflects the choice of the link function and linear predictor using deviance residuals (Nelder & Wedderburn 1972). McCullagh and Nelder (1989) state that through the appropriate link function and linearity of the systematic component, the desired error distribution of the deviance residuals can be achieved which should resemble normal theory residual plots, except for certain plots, like in the case of binomial errors. Standardised deviance residuals are approximately normal which is preferable to Pearson residuals that tend to reflect any skewness of the underlying distribution. Plotting the standardised deviance residuals against the fitted values can provide an informal check of the goodness of fit depending on the type of generalised linear model, where any curvature could suggest the incorrect choice of link function, omitted independent variables or the omission of quadratic terms in the independent variables (Davidson & Snell 1991).

The selection of generalised linear models in this study involved choosing the appropriate distribution of Y and choosing the relationship between η and μ. Three candidate combinations of model families and link functions were fit to the data, specifically the gamma log model, the normal log model and the lognormal identity model. These models are hereafter referred to as the gamma model, the normal model and lognormal model.

Prices of adjacent properties are often related which can lead to correlation in the residuals of regression models, violating the assumption of independence (Bourassa et al. 2007). Spatial autocorrelation or dependency is a challenging problem in real estate modelling where correlation manifests in two-dimensional space unlike serial correlation which is one dimensional. Bourassa et al. (2007) found that the inclusion of a submarket dummy variable accounted for spatial autocorrelation and outperformed geostatistical and lattice approaches. A similar approach was adopted in this study where a factor area variable was included to account for spatial dependence in listing prices.

Variograms were utilised to understand the spatial autocorrelation structure. Variograms display the dissimilarity of observations that vary in space as a function of the distance between them (Ploner 1999). The sill represents spatially autocorrelated sample locations, and the range is where the distance flattens out and the sample locations are no longer spatially autocorrelated. A variogram will be flat when no correlation or low correlation is present which indicates randomness in the structure (Chiles & Delfiner 1999). The nugget effect is an important concept in variograms and describes the variability between observations that are closely spaced which could be inherent in the data or because of the sampling component (Clark 2010). Therefore, in the context of this study, a large nugget effect could be the product of closely clustered properties with similarly signed and order of magnitude residuals that would overestimate the amount of spatial dependency. A prevalent test developed by Moran (1950) is a two-dimensional specification test for spatial autocorrelation, analogous to a test of univariate time series correlation (Anselin 2006): I=e′We/S0(e′ e)/n[Eqn 7] where e represents the regression model residuals and w is the spatial weighting matrix and s₀ is the standardisation factor that relates to the sum of weights for the nonzero cross-products. This test is applied to test for the presence of spatial autocorrelation.

The open source programming language R was used to perform the statistical analysis in this research (R Core Team 2018). The data were provided by an online residential property portal that aggregates listings from real estate agencies throughout South Africa. The period of the data is from January 2013 to August 2017, and Table 1 describes the data set.

The data were inspected for consistency where categorical variables were standardised and data types were constrained to the correct type. Property listings were duplicated by different agencies resulting in many properties being captured more than once. This could result in biased estimates, and therefore duplicate properties were identified and removed using row-wise string matching. Missing values resulted in the removal of observations. The summary statistics show that the spread of the numeric characteristic variables was large with lower and upper bounds that are unlikely.

Real estate agents populate data into automated feeds which could result in incorrect data capturing and anomalous data. To deal with the anomalous data, an autoencoder was developed using the H2O open source machine learning framework (Ledell et al. 2019). An autoencoder is a deep learning neural network aimed at reducing the feature space which can be viewed as a non-linear alternative to principal component analysis (Hastie, Tibshirani & Wainwright 2015). Given enough data, the network will learn the identity of the data via non-linear reduced representation of the original data (Candel et al. 2018). A high reconstruction error for a data point indicates that the data point does not match the learned pattern and is anomalous. Lower limits were set on the listing price, size and lot variables. Listing price was set to ≥ ZAR 200 000; size and lot were set to ≥ 35 m². These figures were chosen based on the ABSA Bank property price index which was used as a guideline (Luus 2002). Properties with a reconstruction mean squared error ≥ 9.39e-07 were discounted. Therefore, based on the results of the autoencoder, properties in the top 5th percentile of the reconstruction error were treated as anomalous.

The final data set consisted of 415 200 properties, and the spread of the variable distributions was greatly reduced after the anomalous data points were removed, evident in Table 2.

Although the lot variable was used to detect anomalous data, it was discounted for modelling purposes as lot was not applicable for the property type apartment and often omitted. Lot was set to the size variable for apartments during the autoencoder learning stage.

The design approach chosen develops separate cross-sectional models for each year in the data which contrasts to building a single model using pooled cross-sectional data. The chosen framework will facilitate the development of a property price index that won’t change previous estimates when future periods are introduced. A pooled period approach would result in new samples being added to the original sample, which would change previous estimations when constructing a residential property price index (De Haan & Diewert 2011). Therefore, for each year, various generalised linear models were developed where listing price was regressed on the physical and locational attributes of residential properties. Tabular and graphical summaries of the model fit for each of the candidate models are presented for comparative purposes. Thereafter, based on generalisability and goodness of fit, the best model is selected and expounded upon.

The data were split into two sets, training and validation, where 70% of the data were used for training and 30% of the data were used for validating the models. This was done to test model generalisability on unseen data for the development of future models. The holdout data for a given model provide a more robust estimate of the generalisation error compared to the training error (Blum, Kalai & Langford 1999). Partitioning the data into training and holdout sets for each year involved writing a function to ensure that the splits were random and that the distribution of the response was similar for each split and to the original data. The function ensured that each area factor level was present in each split. Model performance and generalisation was tested using the root mean squared error (RMSE) which is a measure of spread that compares the closeness of the model outcomes to the observed data (Gujarati 2004). A lower RMSE is indicative of less variability between model estimates and the observed data. The Akaike information criterion (AIC) statistics were also computed. When comparing models, the AIC is useful for model selection as it provides an assessment of the quality of different models given a set of data (Greene 2003). A lower AIC is indicative of better fit. Akaike information criterion concomitantly considers goodness of fit using the likelihood function whilst penalising model complexity through the number of parameters. Model selection was based on a combination of reported statistics, namely deviance explained, holdout RMSE, AIC and model fit based on diagnostic residual plots. Table 3 details the results of each yearly model fit.

Each model produced consistent deviance explained statistics for each year respectively, where the gamma and lognormal models shared the highest amount of deviance explained. Moreover, the gamma and lognormal models appear very similar in terms of holdout RMSE and AIC statistics. The AICs produced by the lognormal models were not directly comparable to the other models as the response variable was on the logarithmic scale. The AICs of the lognormal models were made comparable by subtracting the sum of logarithms of the response variable from the likelihood. Based solely on the AICs, the lognormal models appear to fit the data the best, as they consistently produced the lowest AIC statistics. Considering only the holdout RMSE statistics, the normal model outperformed the two other models with consistently lower RMSE statistics each year. No evidence of overfitting is present as the training and holdout RMSEs are quite similar, indicating that the models generalise to unseen data. This suggests model robustness to the introduction of future periods.

Discerning the best model based solely on the goodness-of-fit measures reported above is difficult, and a graphical examination of the residuals is necessary. The goodness-of-fit residual diagnostic plots for each yearly model are illustrated in Figures 1 and 2 from left to right, beginning with the gamma model, followed by the normal model and finally the lognormal model. Figure 1 presents the residuals versus fitted values where the y-axis represents the deviance residuals and x-axis represents the fitted values. Figure 2 presents the quantile–quantile (Q–Q) plots for normality.

The fitted versus residual diagnostic plots for the gamma and lognormal models are very similar and do not indicate any discernible pattern in the deviance residuals, one of the required assumptions. However, the normal model shows signs of heteroscedasticity at the upper quantiles, violating the assumption of constant variance.

None of the plots is perfectly normal with deviation at the upper and lower quintiles; the S-shaped curves indicate heavy tailed residual distributions. The gamma and lognormal Q–Q plots appear the best behaved in terms of the normality assumption. The normal model appears to fit the data poorly in terms of the diagnostic plots, whilst the gamma and lognormal models appear to represent the data the best.

A possible caveat of using the lognormal model for modelling listing prices is that the expected values are on the log scale and back transformation is necessary. Transforming expected values from the log scale back to the original scale by means of exponentiation results in geometric mean estimates and not arithmetic mean estimates (Olivier et al. 2008). However, the natural logarithm is monotonic, and the back transformed estimates are equivalent to median estimates if the distribution of log(x) is symmetric (Musset 2006). An appealing feature of the gamma and normal models is that expected values are kept on the original scale where arithmetic mean expected values are computed. For this reason, the lognormal models are discounted from the candidate model selection. The gamma models are chosen over the normal models based on the diagnostic plots, lower AICs and similar holdout RMSEs. A discussion of the gamma modelling results ensues.

The property type factor variable included six levels, namely apartment, cluster, duplex, house, simplex and townhouse. The property type apartment was used as the reference level, resulting in the other property types being compared to this level. Table 4 tabulates the beta coefficient estimates for each covariate along with the corresponding p-values. To make reporting succinct, Table 4 discounted the area (factor variable) coefficients as there were over 2000 factor levels present in the data that varied between years. The area variable was used as a control variable to account for variability amongst listing prices and to account for the spatial dependency in the data.

The coefficients given by β^ are expressed as percentage effects. The covariates log(Size) and the number of bathrooms were consistently statistically significant for each year. The natural logarithm was applied to the size covariate to improve linearity where Figure 3 shows no discernible pattern in the plots indicating this transform was appropriate. The coefficients can be interpreted as follows:

It is evident from Table 4 that each additional bathroom, on average, contributes more to the listing prices of homes than each additional bedroom. An appealing feature of this parametric framework is the transparency and interpretability of the model coefficients. Property market participants are able to make informed decisions about renovating their homes or making comparative buying decisions by examining the marginal utility of different characteristics.

Plotting the residuals against individual covariates of the linear predictor should result in a null pattern, like the residual versus fitted values plot (McCullagh & Nelder 1989). The natural logarithm was applied to the size covariate to improve linearity where Figure 4 shows no discernible pattern in the plots indicating this transform was appropriate.

The analysis of deviance presented in Table 5 indicates that the residual deviance for each yearly gamma model was consistently lower than null deviance. This means that the covariates accounted for greater deviance explained than intercept only models and as such indicates a good fit.

The modelling of spatial data in this study required the assessment of the assumption of independence, which was investigated using several plots and by performing hypothesis tests. Variograms quantify the spatial dependence in data by describing the spatial variance. The yearly gamma models’ residuals were plotted using spherical variograms and are presented in Figure 4 where similarities were found between all models. The ranges, distances beyond which the data are no longer correlated, are quite long, which suggests spatial autocorrelation is not an issue in the modelling results. The nugget effects as a percentage of the total sills are quite large which could indicate some variation at a small scale.

A permutation test for Morans I was applied to formally test for the presence of spatial autocorrelation where, under the null hypothesis, the data are randomly dispersed. The Morans I statistic or correlation coefficient ranges between -1 and 1, where -1 shows perfect negative spatial autocorrelation and 1 shows perfect positive spatial autocorrelation. Hundreds of permutations were run, 999 in total, for each yearly gamma model. The results of the tests are presented in Table 6 which indicate a weak negative correlation. Formally, at an alpha of 0.05, there is not enough evidence to reject the null hypothesis of no spatial autocorrelation for each yearly gamma model. This coincides with the findings of Bourassa et al. (2007) where the addition of a location dummy variable accounted for spatial dependence adequately.

Residential property is a barometer of individual and collective wealth and acts as measure of financial stability in an economy. Measuring residential property prices is difficult because of the heterogeneity thereof. The estimation of residential property prices using hedonic modelling is pervasive in real estate economic literature where log linear models are typically employed. This article investigated generalised linear models as an alternative to log linear models to develop hedonic price functions to estimate residential property listing prices in South Africa over a 5-year period. The gamma generalised linear model provided the best fit and good generalisability whilst keeping the expected values on the original scale, which is an appealing alterative to log linear models. The spatial dependence of residential properties was effectively accounted for by including an area factor variable, supported by variograms and Morans I permutation tests, showing no evidence to reject the null hypothesis of no spatial autocorrelation. This framework provides property market participants with the ability to quantify the utility derived over the marginal distribution of the physical characteristics of residential properties. This research presents the groundwork to create a property price index where index number theory could be applied to the hedonic price models to measure price inflation over time.