THE "HEDONIC" MODEL USED IN THE ODPM'S HOUSE PRICE INDICES

A description of the regression model developed for the new ODPM House Price Index that was launched in September 2003

The new monthly House Price Index published by ODPM is based on a multivariate fixed effects regression model ("hedonic" model) of house prices. This paper provides a technical description of the underlying model and how it may be used to estimate average house prices, house price movements and the house price index. For clarity of exposition, the paper contains the following sections:

1. Introduction
2. Specification of the Hedonic Regression Model
3. Estimation of the Expected House Price for a Cell
4. Estimation of the House Price Index
5. Weighted Regressions

Appendix - Estimation of Cumulants

1. Introduction

The rationale for using the hedonic regression approach is that the number of "cells" (that is, possible combinations of different sets of values for the main effects considered) is so large that it is not possible to estimate accurately, or at all, average house prices for each cell individually from the data available. We assume that the influence of each main effect or interaction is the same across all cells and we include in the model those main effects and interactions which "best" explain variation in house prices. The concept of "best" is determined according to Schwarz's Bayesian Criterion (see Schwarz, 1978 and Shi & Tsai, 2002).

The model used includes the following seven main effects:

- location (local authority district or London borough)
- local authority cluster (an ONS classification of local authorities)
- type of neighbourhood (ACORN)
- dwelling type
- number of habitable rooms (or bedrooms)
- old/new
- first-time buyer/former owner occupier (FTB/FOO)

and the following three interactions:

- ACORN × dwelling type
- ACORN × FTB/FOO
- dwelling type × old/new.

For ease of notation, we denote each cell by a single subscript $i$ (or $j$). Associated with each cell $i$ is a vector $x\_{yi}$ containing the values of the variables in the regression model for cell $i$. The subscript $y$ indicates the calendar year to which the data relate. This is because the new house price index is annually chain-linked and we need to allow for the possibility that the model may change at the January link months (because more housing characteristics become available, for example).

Note that, because all the main effects in the model are categorical, they are represented by dummy, dichotomous covariates and the vector $x\_{yi}$ contains only the values 1 or 0, according to whether the relevant characteristic is present or not in the cell. Note also that some elements of $x\_{yi}$ are mutually associated (for example, a dwelling type cannot be both detached and terraced) and all covariates must be linearly independent, to ensure the validity of matrix inversion in estimation. Thus the main effect old/new will be represented by a single covariate with the value (for example) 0 for old and 1 for new. 2. Specification of the Hedonic Regression Model

In month $m$ of year $y$, we assume the following model for the price $P\_{ymik}$ of dwelling $k$ in cell $i$:

$$L\_{ymik} = \\ln(P\_{ymik}) = x\_{yi}'\\beta\_{ym} + \\epsilon\_{ymik} \\quad 2.1$$

where: $\\beta\_{ym}$ is a column vector of $p_y$ unknown parameters which determine the proportionate impact on the expected price of the $p_y$ dichotomous covariates included in the column vector $x\_{yi}$;

and the $\\epsilon\_{ymik}$ are independent, identically distributed random variables with the following properties:

$$E[\\epsilon\_{ymik}] = 0; \\quad \\text{Var}[\\epsilon\_{ymik}] = \\sigma\_{ym}^2; \\quad \\text{Cov}[\\epsilon\_{ymik}, \\epsilon\_{znjl}] = 0 \\quad \\text{for } ymik \\neq znjl \\quad 2.2$$

We also assume that all the higher order central moments of $\\epsilon\_{ymik}$ exist. This is a necessary, but plausible, assumption to cater for the observed circumstance that the $\\epsilon\_{ymik}$ are not Normally distributed. In particular, we assume:

$$E[\\epsilon\_{ymik}^3] = \\mu\_{3ym} = \\kappa\_{3ym}; \\quad E[\\epsilon\_{ymik}^4] = \\mu\_{4ym} = \\kappa\_{4ym} + 3\\sigma\_{ym}^4$$

The expected price of a dwelling in cell $i$, in month $m$ of year $y$ is then:

$$P\_{ymi} = E[P\_{ymik}] = E[\\exp(L\_{ymik})] = \\exp(x\_{yi}'\\beta\_{ym})E[\\exp(\\epsilon\_{ymik})]$$

From the standard theory of moment generating functions (see, for example, chapter 3 of Stuart & Ord, 1987), we have, for random variable $x$ and auxiliary variable $t$:

$$\\ln(E[\\exp{xt}]) = \\sum\_{r=1}^{\\infty} \\frac{\\kappa_r t^r}{r!} \\quad \\text{That is: } E[\\exp{xt}] = \\exp\\left(\\sum\_{r=1}^{\\infty} \\frac{\\kappa_r t^r}{r!}\\right) \\quad 2.3$$

where $\\kappa_r$ is the $r^{th}$ cumulant of the distribution of $x$.

For the house price data and models examined, $\\kappa_2^2 = \\sigma^2 = 0.15$, $\\kappa_3 = 0$, $\\kappa_4 = 0.03$. Higher order cumulants are negligible and may be ignored. For $x = \\varepsilon\_{ymik}$ and $t=1$, we therefore have:

$$E\\left[ \\exp(\\varepsilon\_{ymik}) \\right] = \\exp\\left( 0 + \\frac{\\sigma\_{ym}^2}{2!} + \\frac{\\kappa\_{3ym}}{3!} + \\frac{\\kappa\_{4ym}}{4!} \\right)$$

$$= \\exp\\left( \\frac{\\gamma\_{ym} \\sigma\_{ym}^2}{2} \\right)$$

where $\\gamma\_{ym} = 1 + \\frac{\\kappa\_{3ym}}{3\\sigma\_{ym}^2} + \\frac{\\kappa\_{4ym}}{12\\sigma\_{ym}^2}$ 2.4

is an adjustment factor to correct for non-Normality.

Thus: $P\_{ymi} = E\\left[ P\_{ymik} \\right] = \\exp\\left( x\_{ym}^T \\beta\_{ym} + \\frac{1}{2} \\gamma\_{ym} \\sigma\_{ym}^2 \\right)$ 2.5 3. Estimation of the Expected House Price for a Cell

In practice, the vector $\\beta\_{ym}$ is not known and has to be estimated. For this purpose, we use the following notation:

- $l\_{ym}$ = column vector of all the $n\_{ym}$ elements $l\_{ymik}$ in month $m$ of year $y$;
- $X\_{ym}$ = matrix of housing characteristics for the sample in month $m$ of year $y$; each row of the matrix is the row vector $x'_{yi}$ of housing characteristics for the corresponding element in vector $l_{ym}$;
- $\\varepsilon\_{ym}$ = column vector of all the $n\_{ym}$ random components $\\varepsilon\_{ymik}$ in month $m$ of year $y$.

Equation 2.1 above relates to a single dwelling. The $n\_{ym}$ equations for month $m$ of year $y$ may be summarised in the matrix equation:

$$l\_{ym} = X\_{ym} \\beta\_{ym} + \\varepsilon\_{ym}$$

with:

$$E[\\varepsilon\_{ym}] = 0; \\quad \\text{Var}[\\varepsilon\_{ym}] = \\sigma\_{ym}^2 I \\quad 3.1$$

Using ordinary least squares estimation, the estimator for $\\beta\_{ym}$ is:

$$\\hat{\\beta}_{ym} = (X'_{ym}X\_{ym})^{-1}X'_{ym}l_{ym} \\quad 3.2$$

with:

$$\\text{Var}(\\hat{\\beta}_{ym}) = \\sigma_{ym}^2 (X'_{ym}X_{ym})^{-1} = V\_{ym} \\quad \\text{(for notational convenience)} \\quad 3.3$$

(see, for example, Johnston, 1972).

In practice, the regression model for the House Price Index uses weighted least squares estimation in order to cater for missing values for some covariates. Section 5 below describes the reasons for this, the treatment applied, the impact on the model and amended formulae for $\\hat{\\beta}_{ym}$ and $\\text{Var}(\\hat{\\beta}_{ym})$.

We may therefore consider the following estimator for $P\_{ymi}$:

$$\\tilde{P}_{ymi} = \\exp(x'_{yi}\\hat{\\beta}_{ym} + \\frac{1}{2} \\gamma_{ym} \\sigma\_{ym}^2) \\quad 3.4$$

where the unknown vector $\\beta\_{ym}$ is replaced by the estimator $\\hat{\\beta}\_{ym}$ from equation 3.2.

However, because the vector $\\hat{\\beta}_{ym}$ is random, the expected value of $\\tilde{P}_{ymi}$ is obtained according to equation 2.3 above. We assume that the scalar random variable $x'_{yi}\\hat{\\beta}_{ym}$ is approximately Normal, because the number of parameters in $\\hat{\\beta}\_{ym}$ is very much less than the number of observations on which the estimates are based. With this approximation, only the first two cumulants are non-zero and we have:

\[ E\[\\tilde{P}_{ymi}\] = \\exp\\left(\\mathbf{x}'_yi \\mathbf{\\beta}_{ym} + \\frac{1}{2} \\text{Var}\\left\[\\mathbf{x}'_yi \\mathbf{\\hat{\\beta}}_{ym}\\right\] + \\frac{1}{2} \\gamma_{ym} \\sigma\_{ym}^2\\right) \]

\[ = \\exp\\left(\\mathbf{x}'_yi \\mathbf{\\beta}_{ym} + \\frac{1}{2} \\mathbf{x}'_yi \\mathbf{V}_{ym} \\mathbf{x}_yi + \\frac{1}{2} \\gamma_{ym} \\sigma\_{ym}^2\\right) \\quad 3.5 \]

The estimator (\\tilde{P}\_{ymi}) is therefore biased upwards by the factor (\\exp\\left(\\frac{1}{2} \\mathbf{x}'_yi \\mathbf{V}_{ym} \\mathbf{x}\_yi\\right)). An approximately unbiased estimator is therefore:

\[ \\hat{P}_{ymi} = \\exp\\left(\\mathbf{x}'_yi \\mathbf{\\hat{\\beta}}_{ym} - \\frac{1}{2} \\mathbf{x}'_yi \\mathbf{V}_{ym} \\mathbf{x}_yi + \\frac{1}{2} \\gamma_{ym} \\sigma_{ym}^2\\right) \\quad 3.6 \]

To estimate (\\sigma\_{ym}^2), (\\gamma\_{ym}) and (\\mathbf{V}_{ym}), we use sample estimates of (\\sigma_{ym}^2), (\\kappa\_{ym}) and (\\kappa\_{4ym}). Because these variance and higher order cumulant terms are also estimated, additional bias correction terms should also, in principle, be included. However, these additional bias correction terms are of the order (\\kappa\_{4ym}/n\_{ym}) or less. In practice, (\\kappa\_{4ym}) itself is small and, since (n\_{ym}) is about 20,000, the additional bias correction term is negligible and may be ignored. Usually, (\\hat{\\gamma}\_{ym}) is of the order 1.01.

The appendix to this paper provides formulae, and their derivations, for estimating (\\kappa\_{3ym}) and (\\kappa\_{4ym}). 4. Estimation of the House Price Index

To calculate the house price index, we calculate for each month a weighted average price for the domain of interest, where the domain may be all dwellings covered by the model or a subset of these dwellings, such as dwellings within a region, old or new dwellings, dwellings of a specific type. We use the estimated mean price in each cell from equation 3.6 and a set of weights ( {w\_{Dyi}} ), fixed for each calendar year and

[ \\sum\_{i \\in D} w\_{Dyi} = 1 ]

with ( D ) denotes the set of cells included within the domain of interest. The weighted average price for domain ( D ) in month ( m ) of year ( y ) is then:

\[ \\bar{P}_{Dym} = \\sum_{i \\in D} w\_{Dyi} \\hat{P}_{ymi} = \\sum_{i \\in D} w\_{Dyi} \\exp \\left( x'_y \\hat{\\beta}_{ym} - \\frac{1}{2} x'_y \\hat{\\gamma}_{ym} x_y + \\frac{1}{2} \\hat{\\gamma}_{ym} \\hat{\\sigma}^2_{ym} \\right) \]

where we assume that ( \\hat{\\gamma}_{ym}, \\hat{\\sigma}^2_{ym} ) are sufficiently good estimators for ( \\gamma\_{ym}, \\sigma^2\_{ym} ) that no further bias correction is required (see the discussion after equation 3.6).

The index movement within year ( y ) is ( \\bar{P}_{Dym} / \\bar{P}_{Dy01} ) and the index relative to February 2002, is:

\[ I\_{Dym} = 100 \\frac{\\bar{P}_{D0213}}{\\bar{P}_{D0202}} \\frac{\\bar{P}_{D0313}}{\\bar{P}_{D0301}} \\cdots \\frac{\\bar{P}_{D,y-1,13}}{\\bar{P}_{D,y-1,1}} \\frac{\\bar{P}_{Dym}}{\\bar{P}_{Dy01}} = 100 \\left( \\frac{\\bar{P}_{D0213}}{\\bar{P}_{D0202}} \\prod\_{x=01}^{y-1} \\frac{\\bar{P}_{Dx13}}{\\bar{P}_{Dx01}} \\right) \\frac{\\bar{P}_{Dym}}{\\bar{P}_{Dy01}} \]

where month 13 of year ( y-1 ) = month 1 of year ( y ).

The index may alternatively be expressed as:

\[ I\_{Dym} = 100 \\frac{1}{\\bar{P}_{D0202}} \\frac{\\bar{P}_{D0213}}{\\bar{P}_{D0201}} \\cdots \\frac{\\bar{P}_{D,y-2,13}}{\\bar{P}_{D,y-2,1}} \\frac{\\bar{P}_{D,y-1,13}}{\\bar{P}_{D,y-1,1}} \\frac{\\bar{P}_{Dym}}{\\bar{P}_{Dy01}} = 100 \\left( \\prod_{x=02}^{y-1} \\frac{\\bar{P}_{Dx13}}{\\bar{P}_{Dx+1,01}} \\right) \\frac{\\bar{P}_{Dym}}{\\bar{P}_{D0202}} \]

The ratios ( \\bar{P}_{Dx13} / \\bar{P}_{D,x+1,1} ) are link factors, which ensure continuity of the chain-linked index by correcting for the change, at every January link month, in the weights and model used in the index.

Note that the constraint ( \\sum\_{i \\in D} w\_{Dyi} = 1 ) is only necessary for the calculation of a weighted average price. Equation 4.2 expresses the house price index as a product of ratios of weighted averages. Each ratio contains the same weights in both numerator and denominator. So, for the house price index, only the relative sizes of the weights are important and the constraint ( \\sum\_{i \\in D} w\_{Dyi} = 1 ) is not necessary. 5. Weighted Regressions

Sometimes, the data we receive from lenders are incomplete, with missing values for some covariates. Values may be missing for all values of a covariate (or covariates) from a single lender or only for a few observations. Rather than reject observations with missing values, we have chosen to keep them in the regression analysis with a special value for each covariate to indicate a missing value.

In doing so, we are assuming that, for a given set of known covariates, the missing values are "missing at random". That is, the true values for the missing covariates may be represented as a random realisation from the population distribution of these covariates, for the given set of known covariates. This is equivalent to the more general assumption that covariates not included in the model at all (for example: garage; state of repair of the property) are "missing at random" and may be represented by a random error term. Provided this assumption is valid, there is no bias in the parameters estimated in this way.

However, the assumption of constant variance for the error terms is no longer valid because, clearly, observations with missing values will have greater variance than observations with complete covariate data.

To resolve this problem, we assume that the error variance for observations with missing data is equal to the basic error variance $\\sigma^2\_{ym}$ for observations with complete data multiplied by a factor greater than 1. The non-uniform variances can then be accommodated by applying a weighted regression, using observation weights equal to the reciprocals of estimates of these factors.

For any given combination of missing covariates, we can estimate these weights, which are less than or equal to 1, from regressions based on complete data only. The estimated observation weight, for a given set of missing covariate values, is calculated as the ratio of the mean squared error from the model including all covariates to the mean squared error from the model including only the non-missing covariates.

It may be desirable to re-calculate these weights every month, using the complete data available for that month. However, this would be time-consuming, when the main operational requirement is to validate the data and produce the current month's index values. We believe that these weights will be stable over time and that small errors in these weights would not have any noticeable impact on the estimated indices. However, these observation weights will be re-assessed annually, to ensure that long-term changes to the factors affecting house prices are picked up and to accommodate any changes to the covariates included in the model.

Using weighted regressions, means that the analyses presented in sections 1 to 3 above need to be modified appropriately. There are no modifications to the underlying principles, because weighted regressions can easily be accommodated by appropriate transformations to the dependent variable and covariates. The required modifications and transformations are covered in standard text books (Johnson, 1972, for example). In summary, the following equations from sections 1 to 3 need to be modified as shown, to accommodate the application of weighted regressions.

**Sections 1 and 2**

The original model in section 1 may be re-written as:

\[ L\_{ymik} = \\ln(P\_{ymik}) = x'_{ym} \\beta_{ym} + \\epsilon\_{ymik} \\quad 5.1 \\equiv 2.1a \]

where the modified error terms have the properties:

\[ E\[\\epsilon\_{ymik}^*\] = 0; \\quad \\text{Var}\[\\epsilon\_{ymik}^*\] = \\frac{\\sigma\_{ym}^2}{\\lambda\_{ymik}}; \\quad \\text{Cov}\[\\epsilon\_{ymik}^*, \\epsilon\_{znjl}^*\] = 0 \\quad \\text{for } ymik \\neq znjl \\quad 5.2 \\equiv 2.2a \]

with ( \\lambda\_{ymik} ) being the appropriate weight for observation ( ymik ).

It is possible to transform this model into a model whose error terms have constant variance by transforming the dependent variable and covariates, as follows:

\[ \\sqrt{\\lambda\_{ymik}} L\_{ymik} = \\sqrt{\\lambda\_{ymik}} \\ln(P\_{ymik}) = \\sqrt{\\lambda\_{ymik}} x'_{ym} \\beta_{ym} + \\epsilon\_{ymik} \\quad 5.3 \\equiv 2.1b \]

so that the error term ( \\epsilon\_{ymik} = \\sqrt{\\lambda\_{ymik}} \\epsilon\_{ymik}^\* ) has the homoscedastic property of equations 2.2.

**Section 3**

\[ E\[\\epsilon\_{ym}^*\] = 0; \\quad \\text{Var}\[\\epsilon\_{ym}^*\] = \\sigma\_{ym}^2 \\Omega\_{ym} = \\sigma\_{ym}^2 W\_{ym}^{-1} \\quad 5.4 \\equiv 3.1a \]

where ( W\_{ym} ) is a diagonal matrix of the weights ( {\\lambda\_{ymik}} ).

Using weighted least squares estimation, the estimator for ( \\beta\_{ym} ) is:

\[ \\hat{\\beta}_{ym} = \\left( X'_{ym} W\_{ym} X\_{ym} \\right)^{-1} X'_{ym} W_{ym} l\_{ym} \\quad 5.5 \\equiv 3.2a \]

\[ \\text{Var}(\\hat{\\beta}_{ym}) = \\sigma^2 \\left( X'_{ym} W\_{ym} X\_{ym} \\right)^{-1} = V\_{ym} \\quad \\text{(for notational convenience)} \\quad 5.6 \\equiv 3.3a \]

(see, for example, Johnston, 1972, chapter 7). References Cook, J, and Weisberg, S. (1982) *Residuals and Influence in Regression*, Chapman & Hall, London Johnston, J. (1972) *Econometric Methods*, 2nd ed., McGraw-Hill, Tokyo Schwarz, G. (1978), Estimating the dimension of a model, *Ann. Statist.* , 6, 461-464. Shi, P., Tsai, C-L. (2002), Regression model selection - a residual likelihood approach, *J. R. Statist. Soc. B*, 64, Part 2, 237-252 Stuart, A. and Ord, J.K. (1987) *Kendall's Advanced Theory of Statistics*, 5th ed., Vol.1, Griffin, London Appendix

Estimation of Cumulants

1. Introduction

The model described in section 1 assumes that the error terms are independent and identically distributed but not that they are Normally distributed. As discussed in section 2, the higher order moments, in the form of cumulants, are relevant for the unbiased estimation of average house prices. This appendix develops estimators for the third and fourth cumulants.

For clarity of exposition, the suffices relating to year, month and cell will be omitted. Suffices will define only the observation number, out of ( n ) observations using a model fitted to ( p ) parameters. For the House Price Index, ( n \\approx 20,000 ) and ( p \\approx 200 ).

The weighted residual from the model for observation ( i ) is defined as:

[ e_i = \\sqrt{\\lambda_i} \\hat{r}\_i = \\sqrt{\\lambda_i} (l_i - x_i' \\hat{\\beta}) = \\sqrt{\\lambda_i} l_i - \\sqrt{\\lambda_i} x_i' \\hat{\\beta} ]

Standard regression theory states that the vector of residuals ( e ) may be expressed as a linear function of the vector of independent, identically distributed error terms ( \\varepsilon ) (see, for example, Johnston, 1972): ( e = (I - H)\\varepsilon = M\\varepsilon ), where ( I ) is the identity matrix and ( H ) is the so-called "hat" matrix ( H = X(X'X)^{-1}X' ).

In particular:

[ e_i = \\sum\_{j=1}^{n} M\_{ij} e_j = \\sum\_{j=1}^{n} (I\_{ij} - H\_{ij}) e_j ]

Note that ( l\_{ij} = 1 ), ( l\_{ij} = 0 ) for ( i \\neq j ) and ( (l\_{ij})' = l\_{ij} ) for all ( i, j, r ).

Also:

[ H\_{ii} \\leq 1, \\quad tr(H) = \\sum\_{i=1}^{n} H\_{ii} = p \\quad \\text{and} \\quad \\sum\_{j=1}^{n} H\_{ij}^2 = H\_{ii} ]

because ( H ) is a symmetric, idempotent matrix of rank ( p ) (see, for example, Cook & Weisberg, 1982). 2. Estimating the Error Variance $\\sigma^2$

We illustrate our approach to developing estimators for cumulants by applying it to produce the standard unbiased estimator for the error variance $\\sigma^2$. We have:

$$E \\left[ \\sum\_{i=1}^{n} e_i^2 \\right] = E \\left[ \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} M\_{ij} e_j \\sum\_{k=1}^{n} M\_{ik} e_k \\right]$$

$$= \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} \\sum\_{k=1}^{n} M\_{ij} M\_{ik} E \\left[ e_j e_k \\right]$$

$$= \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} M\_{ij}^2 \\sigma^2$$

$$= \\sigma^2 \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} \\left( I\_{ij} - H\_{ij} \\right)^2$$

$$= \\sigma^2 \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} \\left( I\_{ij} - 2I\_{ij} H\_{ij} + H\_{ij}^2 \\right)$$

$$= \\sigma^2 \\sum\_{i=1}^{n} \\left( 1 - 2H\_{ii} + H\_{ii} \\right)$$

$$= \\sigma^2 (n - p)$$

Hence we have the standard, unbiased estimator for the error variance:

$$\\hat{\\sigma}^2 = \\frac{\\sum\_{i=1}^{n} e_i^2}{(n - p)}$$

3. Estimating the Third Cumulant $\\kappa_3$

$$E \\left[ \\sum\_{i=1}^{n} e_i^3 \\right] = E \\left[ \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} M\_{ij} e_j \\sum\_{k=1}^{n} M\_{ik} e_k \\sum\_{l=1}^{n} M\_{il} e_l \\right]$$

$$= \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} \\sum\_{k=1}^{n} \\sum\_{l=1}^{n} M\_{ij} M\_{ik} M\_{il} E \\left[ e_j e_k e_l \\right]$$

$$= \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} M\_{ij}^3 \\mu_3$$

because the independence of the $e_i$ means that $E[e_j e_k e_l]=0$ if any subscript $j$, $k$ or $l$ is different from the other two. Hence, noting also that $\\mu_3 = \\kappa_3$:

$$E \\left[ \\sum\_{i=1}^{n} e_i^3 \\right] = \\kappa_3 \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} (I\_{ij} - H\_{ij})^3$$

$$= \\kappa_3 \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} \\left( I\_{ij} - 3I\_{ij}H\_{ij} + 3I\_{ij}H\_{ij}^2 - H\_{ij}^3 \\right)$$

$$= \\kappa_3 \\left{ \\sum\_{i=1}^{n} \\left( 1 - 3H\_{ii} + 3H\_{ii}^2 \\right) - \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} H\_{ij}^3 \\right}$$

$$= \\kappa_3 \\left{ (n - 3p) + 3\\sum\_{i=1}^{n} H\_{ii}^2 - \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} H\_{ij}^3 \\right}$$

The mean value of $H\_{ii}$ is $p/n$. If the absolute deviations of the $H\_{ii}$ from the mean $p/n$ are mainly small (that is, there are not too many large values of $H\_{ii}$), then each $H\_{ii}$ will be of the order $p/n$ and $\\frac{p^2}{n} \\ll \\sum\_{i=1}^{n} H\_{ii}^2 \\leq p$. In practice, for the house price data, $\\sum\_{i=1}^{n} H\_{ii}^2 = 4$. This is very small relative to $(n-3p)$ and may be neglected.

Using the same reasoning, the mean value of $H\_{ij}^2$ ($i \\neq j$) is approximately $\\frac{p - p^2/n}{n^2 - n} \\approx \\frac{p}{n^2}$. Again, if the deviations of $|H\_{ij}|$ from the mean absolute value $(\\sqrt{p})/n$ are sufficiently small, we would expect $\\sum\_{i=1}^{n} \\sum\_{j=1}^{n} |H\_{ij}| \\approx \\frac{p^{1.5}}{n}$. In fact, some of the $H\_{ij}$ may be negative and $\\sum\_{i=1}^{n} \\sum\_{j=1}^{n} H\_{ij}^3$ will be smaller than this. So this term, too, may be neglected.

We therefore have the following, approximately unbiased estimator for the third cumulant:

$$\\hat{\\kappa}_3 = \\frac{\\sum_{i=1}^{n} e_i^3}{n - 3p}$$ 4. Estimating the Fourth Cumulant $\\kappa_4$

$$E \\left[ \\sum\_{i=1}^{n} e_i^4 \\right] = E \\left[ \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} M\_{ij} e_j \\sum\_{k=1}^{n} M\_{ik} e_k \\sum\_{l=1}^{n} M\_{il} e_l \\sum\_{m=1}^{n} M\_{im} e_m \\right]$$

$$= \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} \\sum\_{k=1}^{n} \\sum\_{l=1}^{n} \\sum\_{m=1}^{n} M\_{ij} M\_{ik} M\_{il} M\_{im} E \\left[ e_j e_k e_l e_m \\right]$$

The independence of the $e_i$ means that $E[e_j e_k e_l e_m] = 0$ if any subscript $j$, $k$, $l$ or $m$ is different from all the other three. So $E[e_j e_k e_l e_m]$ is not zero only if pairs of subscripts are equal, that is:

$$j=k \\text{ and } l=m \\ (j \\neq l); \\ j=l \\text{ and } k=m \\ (j \\neq k); \\ j=m \\text{ and } k=l \\ (j \\neq k); \\text{ or } j=k=l=m.$$

The first three conditions are equivalent, corresponding merely to rotation of subscripts, so:

$$E \\left[ \\sum\_{i=1}^{n} e_i^4 \\right] = \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} M\_{ij}^4 E \\left[ e_j^4 \\right] + 3 \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} \\sum\_{k=1}^{n} M\_{ij}^2 M\_{ik}^2 E \\left[ e_j^2 e_k^2 \\right]$$

$$= \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} M\_{ij}^4 \\mu_4 + 3 \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} \\sum\_{k=1}^{n} M\_{ij}^2 M\_{ik}^2 \\sigma^4 \\quad \\text{(from the independence of } e_j \\text{ and } e_k \\text{)}$$

$$= \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} M\_{ij}^4 \\kappa_4 + 3 \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} \\sum\_{k=1}^{n} M\_{ij}^2 M\_{ik}^2 \\sigma^4$$

$$= \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} M\_{ij}^4 \\kappa_4 + 3 \\left( \\sum\_{j=1}^{n} M\_{ij}^2 \\right)^2 \\sigma^4$$ Hence:

\[ E \\left[ \\sum\_{i=1}^{n} e_i^4 \\right] = \\kappa_4 \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} (I\_{ij} - H\_{ij})^4 + 3\\sigma^4 \\sum\_{i=1}^{n} \\left{ \\sum\_{j=1}^{n} (I\_{ij} - H\_{ij})^2 \\right}^2 \]

[ = \\kappa_4 \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} \\left( I\_{ij} - 4I\_{ij}H\_{ij} + 6I\_{ij}H\_{ij}^2 - 4I\_{ij}H\_{ij}^3 + H\_{ij}^4 \\right) + 3\\sigma^4 \\sum\_{i=1}^{n} \\left{ \\sum\_{j=1}^{n} (I\_{ij} - 2I\_{ij}H\_{ij} + H\_{ij}^2) \\right}^2 ]

[ = \\kappa_4 \\left{ \\sum\_{i=1}^{n} \\left( 1 - 4H\_{ii} + 6H\_{ii}^2 - 4H\_{ii}^3 \\right) + \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} H\_{ij}^4 \\right} + 3\\sigma^4 \\sum\_{i=1}^{n} \\left{ 1 - H\_{ii} \\right}^2 ]

[ = \\kappa_4 \\left{ n - 4p + 6\\sum\_{i=1}^{n} H\_{ii}^2 - 4\\sum\_{i=1}^{n} H\_{ii}^3 + \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} H\_{ij}^4 \\right} + 3\\sigma^4 \\left{ 1 - 2H\_{ii} + H\_{ii}^2 \\right} ]

[ = \\kappa_4 \\left{ n - 4p + 6\\sum\_{i=1}^{n} H\_{ii}^2 - 4\\sum\_{i=1}^{n} H\_{ii}^3 + \\sum\_{i=1}^{n} \\sum\_{j=1}^{n} H\_{ij}^4 \\right} + 3\\sigma^4 \\left( n - 2p + \\sum\_{i=1}^{n} H\_{ii}^2 \\right) ]

Using the same arguments presented above for the third cumulant ( \\kappa_3 ), we may ignore the summation terms in the line above to obtain the approximation:

\[ E \\left[ \\sum\_{i=1}^{n} e_i^4 \\right] \\approx \\kappa_4 \\left{ n - 4p \\right} + 3\\sigma^4 \\left( n - 2p \\right) \]

from which we have the following, approximately unbiased estimator for the fourth cumulant ( \\kappa_4 ):

\[ \\hat{\\kappa}_4 = \\frac{\\sum_{i=1}^{n} e_i^4 - 3\\hat{\\sigma}^4 (n - 2p)}{(n - 4p)} \]

Note that ( \\hat{\\sigma}^4 = (\\hat{\\sigma}^2)^2 ) is a biased estimator of ( \\sigma^4 ) but the bias is of the order ( \\kappa_4/n ), which is negligible to the degree of approximation used here.