%%%%%%%%%%%%%%%%%%%% author.tex %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% sample root file for your "contribution" to a contributed volume
%
% Use this file as a template for your own input.
%
%%%%%%%%%%%%%%%% Springer %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%% RECOMMENDED %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass[graybox]{svmult}

% choose options for [] as required from the list
% in the Reference Guide

\usepackage{type1cm}        % activate if the above 3 fonts are
                            % not available on your system
%
\usepackage{makeidx}         % allows index generation
\usepackage{graphicx}        % standard LaTeX graphics tool
                             % when including figure files
\usepackage{multicol}        % used for the two-column index
\usepackage{multirow}
\renewcommand{\arraystretch}{1.5}
\usepackage[bottom]{footmisc}% places footnotes at page bottom

\usepackage[utf8x]{inputenc}

\usepackage{caption}
\captionsetup{justification=justified}
\usepackage{tablefootnote}
\usepackage[hidelinks]{hyperref}
\usepackage{url}
\usepackage{marvosym}

\usepackage{apacite}

\AtBeginDocument{\urlstyle{APACsame}}  % so that URLs are same font not monochrome
\bibliographystyle{apacite}

\usepackage{newtxtext}       % 
\usepackage{newtxmath}       % selects Times Roman as basic font

% added commands and userpackages needed for this chapter
\newcommand{\plim}{plim}
\newcommand{\tr}{tr}
\newcommand{\diag}{diag}
\newcommand{\rank}{rank}
\usepackage{eurosym}
%\usepackage{natbib}



% see the list of further useful packages
% in the Reference Guide

\makeindex             % used for the subject index
                       % please use the style svind.ist with
                       % your makeindex program

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\setcounter{chapter}{11} %Set counter to one less than your chapter number.
%%%%%%%%%%%%%%%%%%%
\begin{document}

\title{Modelling Housing Using Multi-dimensional Panel Data}
% Use \titlerunning{Short Title} for an abbreviated version of
% your contribution title if the original one is too long
\author{Badi H. Baltagi and Georges Bresson}
% Use \authorrunning{Short Title} for an abbreviated version of
% your contributing authors if the original one is too long
\authorrunning{Baltagi and Bresson}
%%If only one or two authors: \authorrunning{Surname Author 1 and surname author 2}

\institute{Badi H. Baltagi \Letter \at Department of Economics and Center for Policy Research, Syracuse University, Syracuse, New York, USA, \email{bbaltagi@maxwell.syr.edu} 
\and
Georges Bresson \at Universit\'{e} Paris Panth\'{e}on-Assas, Paris, France, \email{georges.bresson@u-paris2.fr}
}




%
% Use the package "url.sty" to avoid
% problems with special characters
% used in your e-mail or web address
%
\maketitle

\abstract*{This chapter surveys housing models using
multi-dimensional panels. While there is a vast literature on housing models
using two-dimensional panel data, there are only few papers using
multi-dimensional panels. This chapter focuses on housing models, residential
mobility and location choice models derived from discrete choice theory
utilizing multi-dimensional panels. Examples  include nested or hierarchical error
components models where a house is located in a street, within a block, within a city, within a county, etc.
This chapter introduces some basic concepts of utility
functions and discrete choice models used for the hedonic functions and the
residential mobility and location choices. Then it surveys some papers on
multi-dimensional models of housing hedonic price functions focusing on
their estimation methods and their main results. This is followed by a
survey of some papers on multi-dimensional models of residential
mobility and location choice as well as surveying a few
papers on dynamic housing models. It shows that both spatial and temporal dimensions in dynamic systems should be included for hedonic housing models and discrete models of residential location in a multi-dimensional framework. But the inclusion of these multiple dimensions greatly complicates the specification and modeling of such systems. Last, the paper concludes with variational Bayesian approximations which are promising future pathways to potentially overcome many problems in applied modelling of housing and illustrate it using hedonic housing estimation for the city of Paris.}


\abstract{This chapter surveys housing models using
multi-dimensional panels. While there is a vast literature on housing models
using two-dimensional panel data, there are only few papers using
multi-dimensional panels. This chapter focuses on housing models, residential
mobility and location choice models derived from discrete choice theory
utilizing multi-dimensional panels. Examples  include nested or hierarchical error
components models where a house is located in a street, within a block, within a city, within a county, etc.
This chapter introduces some basic concepts of utility
functions and discrete choice models used for the hedonic functions and the
residential mobility and location choices. Then it surveys some papers on
multi-dimensional models of housing hedonic price functions focusing on
their estimation methods and their main results. This is followed by a
survey of some papers on multi-dimensional models of residential
mobility and location choice as well as surveying a few
papers on dynamic housing models. It shows that both spatial and temporal dimensions in dynamic systems should be included for hedonic housing models and discrete models of residential location in a multi-dimensional framework. But the inclusion of these multiple dimensions greatly complicates the specification and modeling of such systems. Last, the paper concludes with variational Bayesian approximations which are promising future pathways to potentially overcome many problems in applied modelling of housing and illustrate it using hedonic housing estimation for the city of Paris.}


\section{Introduction}

\label{CH12X_section1}

This chapter surveys \index{housing!model}housing models using \index{multi-dimensional panels} multi-dimensional
panels. For more than a decade, a huge literature within the
\index{New Economic Geography}New Economic Geography has emerged to study the causes of temporal and spatial variations in \index{housing!price}house prices,
\index{residential!mobility}residential mobility and \index{location choice}location choice. These are major
household decisions connected with many activities and travel aspects of
households lives. These concepts have been widely researched in various
fields including economics, sociology, geography, urban planning,
transportation, etc. \index{location choice}Location choices and \index{housing!investments}housing investments are inherently
dynamic decisions. Moreover, the choice for a household to locate in a given
area is a complex decision that is influenced by, among other things, the
structural elements of a \index{dwelling}dwelling as well as the property's spatial
relationship to certain \index{amenities}amenities. One source of \index{spatial!heterogeneity}spatial heterogeneity comes from the natural hierarchical and
\index{nested!structure}nested structure of the locations of houses, located in a street, within a
block, within a city, within a county, within a region, etc. There is a vast
literature on such topics mainly using time series and longitudinal
(two-dimensional (2D)) data but only few papers using a multi-dimensional
(\index{three-dimensional}three-dimensional (3D) and more) framework. In this chapter, we will focus
on \index{housing!model}housing models, \index{residential!mobility}residential mobility and \index{location choice}location choice models derived
from \index{discrete choice}discrete choice theory focusing on examples that use multi-dimensional
panels.

For example, \citeA{Baltagi2014} focus on the estimation of \index{housing!model example!UK}UK \index{housing!price}house prices
in which \index{spatio-temporal}spatio-temporal variations in \index{housing!price}house prices are driven by supply and
demand conditions, with \index{spatial!effect}spatial effects coming from two distinct sources.
One is the direct dependence of \index{housing!price}house prices in a given locality on house
prices in nearby localities. The second source of \index{spatial!heterogeneity}spatial heterogeneity
comes from the presence of hierarchical error components which represent the
impact of local (district) effects embedded within wider (county) effects.
The panel data includes 353 local authority districts in\ England over the
period 2000-2007. This is done using \index{instrumental variable}instrumental variable
estimation. Another example is \citeA{Baltagi2015} who estimate a \index{hedonic!housing model}hedonic housing model based on flats sold in the city of \index{housing!model example!Paris}Paris
over the period 1990-2003. This is done using \index{maximum likelihood}maximum likelihood
estimation, taking into account the \index{nested!structure}nested structure of the data. \index{housing!model example!Paris}Paris is
historically divided into 20 \textit{arrondissements}, each divided into
four \textit{quartiers} (quarters), which in turn contain between 15 and 169
blocks (\textit{\^{\i}lot}, in French) per \textit{quartier}.

In Sect. \ref{CH12X_section2} we introduce some basic concepts of utility
functions and \index{discrete choice}discrete choice models used for the \index{hedonic!housing model}hedonic functions, the
\index{residential!mobility}residential mobility and \index{location choice}location choices. Section \ref{CH12X_section3} deals
with multi-dimensional models of housing \index{hedonic!price}hedonic price functions, their
estimation methods and some results. Section \ref{CH12X_section4} analyses
some multi-dimensional models of \index{residential!mobility}residential mobility and \index{location choice}location choice.
Section \ref{CH12X_section5} focuses on multi-dimensional \index{dynamic model}dynamic models of
\index{housing!model}housing models. Section \ref{CH12X_section6} highlights \index{variational inference}variational Bayesian inference, and more specifically \index{variational inference!mean field variational Bayes approximation (MFVB)}mean field variational Bayes approximations to specify and estimate complex multi-dimensional housing models while Section \ref{CH12X_section7} concludes.

\section{Discrete Choice Models and Hedonic Price Functions: A Quick Overview%
}

\label{CH12X_section2}

The pioneering work by Daniel McFadden on \index{location choice}location choice is an obvious
starting point for a discussion on \index{housing!model}housing models. One generally considers a
household $i$ who chooses to locate in \index{neighborhood}neighborhood $j$ and buy house type $k
$. A standard \index{random!utility}random utility model (see \citeA{Holmes2014} for
instance) assumes that the \index{indirect utility}indirect utility of household $i$ for
location $j$ and house $k$ is given by
\begin{equation*}
u_{ijk}=X_{j}^{\prime }\beta +Z_{k}^{\prime }\gamma +(y_{i}-p_{jk})\alpha
+\varepsilon _{ijk}=f_{ijk}(.)+\varepsilon _{ijk}\,,  \label{CH12X_eq_1}
\end{equation*}%
where $X_{j}$ is a vector of observed characteristics of location $j$, $Z_{k}
$ is a vector of observed characteristics for house $k$, $y_{i}$ is the
household income and $p_{jk}$ is the price of housing type $k$ in location $j
$. Each household chooses the \index{neighborhood}neighborhood-housing pair that
maximizes utility. Under the assumption that the error terms $\varepsilon
_{ijk}$ are independent and identically distributed ($i.i.d.$) across $i$, $j
$ and $k$ and follow a type I extreme value distribution, \citeA{McFadden1973}
(see also \citeA{McFadden1974, McFadden1978}), derived the well-known \index{conditional logit}conditional logit choice probabilities:
\begin{equation*}
Pr\left[ d_{ijk}=1\right] =\frac{\exp \left( f_{ijk}(.)\right) }{%
\sum_{j=1}^{J}\sum_{k=1}^{K}\exp \left( f_{ijk}(.)\right) }\,,
\label{CH12X_eq_2}
\end{equation*}%
where $d_{ijk}=1$ if household $i$ has chosen \index{neighborhood}neighborhood $j$ and house
type $k$ and zero otherwise. But, the \index{independence of irrelevant alternatives}independence of irrelevant alternatives (IIA) property of this model is unattractive. \citeA
{McFadden1978} proposed the use of a \index{generalized!extreme value distribution}generalized extreme value distribution
for the error terms which gives rise to the \index{nested!logit}nested logit model and allows
one to relax the assumption that idiosyncratic tastes are independent across
locations and houses. But, we need to choose the \index{nested!structure}nesting structure before
estimation, mainly if the \index{nested!structure}nested structure is not natural and if we do not
have knowledge about the \index{neighborhood}neighborhood structure. One solution is to use
\index{random!coefficients}random coefficients $\beta _{i}$, $\gamma _{i}$ and $\alpha _{i}$ instead of
fixed coefficients $\beta $, $\gamma $ and $\alpha $. Estimation with \index{random!coefficients}random coefficients is challenging and needs the use of \index{simulation-based estimator}simulation-based estimators
(SBE) (see \citeA{Newey1974} or \citeA{Judd1998}).

Moreover, \index{Bayesian!estimators}Bayesian estimators are also well suited for the
estimation of \index{discrete choice}discrete choice models with \index{random!coefficients}random coefficients.
One application of such model with SBE has been done by \citeA
{Hastings2006} who study the effects of open enrollment policies under a
particular parent choice mechanism, \index{sorting}sorting households among schools within
the Mecklenburg Charlotte school district, \index{housing!model example!US North Carolina}North Carolina. \citeA{Bajari2005}
used \index{Bayesian!methods}Bayesian methods to study \index{housing!demand}housing demand explaining racial
segregation in cities.

Demand estimation has also focused on the role of unobserved
\index{neighborhood}neighborhood characteristics or \index{housing!quality}housing quality $\zeta _{j}$. In that case,
the indirect utility function is written as
\begin{equation*}
u_{ijk}=X_{j}^{\prime }\beta +Z_{k}^{\prime }\gamma +(y_{i}-p_{jk})\alpha
+\zeta _{j}+\varepsilon _{ijk}\,.  \label{CH12X_eq_3}
\end{equation*}%
Unobserved \index{neighborhood}neighborhood characteristics can be recovered by matching the
observed market shares of community $j$. Then, the remaining parameters can
be estimated by a \index{generalized!method of moments}generalized method of moments (GMM) estimator
using \index{instrumental variable}instrumental variables (IV) to deal with the correlation between
\index{housing!price}housing price $p_{jk}$ and unobserved \index{neighborhood}neighborhood characteristics or
housing \index{quality}quality $\zeta _{j}$. \citeA{Bayer2007}, using a two-dimensional (2D)
panel data, estimate household preferences for school and \index{neighborhood}neighborhood
\index{attributes}attributes in the presence of \index{sorting}sorting. The model embeds a boundary
discontinuity design in a heterogeneous \index{residential!choice}residential choice model, addressing
the \index{endogeneity}endogeneity of the school and \index{neighborhood}neighborhood characteristics. Their
application concerns a restricted-access version of the 1990 \index{U.S. Census}U.S. Census, that links detailed characteristics for nearly a quarter of a
million households and their houses in the \index{housing!model example!US San Francisco}San Francisco Bay Area with their
precise \index{residential!location}residential location. \citeA{Bayer2016}, using a \index{three-dimensional}three-dimensional
panel data (3D), develop a \index{dynamic model}dynamic model of \index{neighborhood}neighborhood choice (see Sect. %
\ref{CH12X_section5}). They capture observed and unobserved preference
\index{heterogeneity}heterogeneity across households and locations of housing transactions in the
\index{housing!model example!US San Francisco}San Francisco Bay Area from 1994 to 2004.

We turn now to \index{hedonic!measure}hedonic measures with a strong theoretical grounding (see,
for example, \citeA{Griliches1971}, \citeA{Rosen1974}, \citeA{Nelson1977}, \citeA
{Blomquist1981}, \citeA{Blomquist1982} among others). Also, the use of
regression techniques to control for compositional and \index{quality!change}quality change (\citeA{Witte1979}, \citeA{Brown1982}, \citeA{Meese1997}, to mention a few).
The \index{hedonic!price}hedonic pricing method is based on the fact that prices of goods (in our
case, houses) in a market are affected by their characteristics. This method
estimates the value of a commodity based on people's \index{willingness!to pay} willingness to pay for the commodity as and when its characteristics change. 
In \index{real estate}real estate economics, hedonic pricing is used to adjust for the problems
associated with looking for a \index{dwelling}dwelling that is as heterogeneous as
buildings. The \index{hedonic!price}hedonic pricing function, which explains the price of a
house, will be affected by among other things, the structural
characteristics of the house, the characteristics of the \index{neighborhood}neighborhood and
the \index{environmental characteristics}environmental characteristics.

Since the seminal work of \citeA{Rosen1974}, one generally uses a two-stage
procedure for estimating the \index{hedonic!price}hedonic price function of the \index{dwelling}dwelling and for
the recovery of marginal \index{willingness!to pay}willingness to pay functions of heterogeneous
individuals for the characteristics of differentiated products. Basically,
hedonic models of \index{housing!price}housing price relate the price (or the logarithm of the
price per square meter) to among other things, the characteristics of the
\index{dwelling}dwellings $p_{jk}=f(Z_{k}^{\prime },...)$. The price gradient associated
with this \index{hedonic!price}hedonic price function $\partial p_{jk}/\partial Z_{kl}$ denotes
the implicit price of the amenity $Z_{kl}$ (number of rooms, \index{quality}quality of air,
etc.). The second stage of Rosen's procedure seeks to recover the
coefficients of demand (or marginal \index{willingness!to pay}willingness to pay) and supply (or
marginal \index{willingness!to accept}willingness to accept) functions for the \index{attributes}attribute $Z_{kl}$ from
the first-order conditions of the equilibrium relationships: $\partial
p_{jk}/\partial Z_{kl}=f_{d}\left( Z_{k},B_{k}\right) $ for demand and $%
\partial p_{jk}/\partial Z_{kl}=f_{s}\left( Z_{k},S_{k}\right) $ for supply
where $B_{k}$ and $S_{k}$ represent \index{attributes}attributes of the buyer and seller of
house $k$. \citeA{Bartik1987} and \citeA{Epple1987} have described a source of
\index{endogeneity}endogeneity in the second stage of Rosen's procedure that is difficult to
overcome without exclusion restriction arguments or the use of IV methods.
This has led researchers to avoid the estimation of marginal
\index{willingness!to pay}willingness to pay functions altogether, relying instead on the first-stage
\index{hedonic!price}hedonic price function and limiting the analysis to the evaluation of
marginal changes in \index{amenities}amenities (see \citeA{Gayer2000}, \citeA{Bishop2011} to
mention a few).

In some studies, \index{dwelling}dwellings were assumed to be stratified into blocks or
communities $j$ where prices are homogeneous and price trends are roughly
parallel. Ideally a model could be estimated in each \index{neighborhood}neighborhood and the
elementary geographic zones could be very small sub-markets. In that case,
each model is estimated in a particular block, all variables are \textit{de
facto} interacted with the block. So, spatial location is not without
consequences and hedonic \index{housing!price}housing price models should incorporate spatial
effects. In the econometric literature, \index{spatial!effect}spatial effects may result from
\index{spatial!dependence}spatial dependence or from \index{spatial!heterogeneity}spatial heterogeneity. \index{spatial!dependence}Spatial dependence means
that observations at location $j$ depend on other observations at locations $%
l\neq j$. \index{spatial!heterogeneity}Spatial heterogeneity refers to variation in relationships over
space and more precisely over every point in space. The distinction is
coming from the structure of the dependence which can be related to location
and distance, both in a geographic space as well as in a more general
economic or social network space (see \citeA{Anselin2001}, \citeA{Anselin2008}%
).

For \index{spatial!effect}spatial effects in \index{real estate}real estate, a lot of housing models have been
estimated in a 2D framework on panel data with two indexes $j$ and $t$
generally for location and time associated with \index{spatial!weight}spatial weight matrices (see
for instance \citeA{Baltagi2011}, \citeA{Bresson2011}, \citeA{Fingleton2008},
\citeA{Glaeser2008}, \citeA{Holly2010} to mention a few). But very few models
have been developed in a \index{three-dimensional}three-dimensional, \index{four-dimensional}four-dimensional, or more in a
panel data setting. In the next section, we present some of these models and
their associated results for these multi-dimensional frameworks.

\section{Multi-dimensional Models of Housing Hedonic Price Functions: Some
Examples}

\label{CH12X_section3}

\citeA{Baltagi2015} estimate a hedonic housing model based on flats sold in
the city of \index{housing!model example!Paris}Paris over the period 1990-2003. This is done using \index{maximum likelihood}maximum likelihood estimation, taking into account the \index{nested!structure}nested structure of
the data. \index{housing!model example!Paris}Paris is historically divided into 20 \textit{arrondissements},
each divided into four \textit{quartiers} (quarters), which in turn contain
between 15 and 169 blocks (\textit{\^{\i}lot}, in French) per \textit{%
quartier}. The data set used is an \index{unbalanced pseudo-panel data}unbalanced pseudo-panel data
containing 156,896 transactions. The \index{real estate}real estate literature emphasizes the
importance of \index{neighborhood}neighborhoods in determining the value of a house or a flat.
While one can try and include as many as possible of the \index{neighborhood}neighborhood
characteristics in the regression to capture these effects, most attempts
may fall short as many \index{neighborhood}neighborhood characteristics are not observed, as in
our case. One simple method of capturing the effect of \index{neighbor}neighbors' prices
used by \citeA{Baltagi2015} is to estimate a \index{spatial!lag}spatial lag regression
equation with \index{time-varying coefficients}time-varying coefficients:
\begin{equation}
p_{taqif}=\lambda _{t}\tilde{p}_{taqif}+Z_{taqif}\beta +\varepsilon _{taqif}%
\text{ , }| \lambda _{t}| <1\,,  \label{CH12X_eq_4}
\end{equation}%
where $t=1,...,T$ for years, $a=1,...,N$ for \textit{arrondissements}, $%
q=1,...,Q_{ta}$ for \textit{quartiers}, $i=1,...,M_{taq}$ for \textit{\^{\i}%
lots} and $f=1,...,F_{taqi}$ for flats. $p$ is the transaction price (in
logs) for flat $f$, in \textit{\^{\i}lot} $i$ \index{nested!structure}nested in \textit{quartier} $q$%
, which in turn is \index{nested!structure}nested in \textit{arrondissement} $a$ at time $t$. $%
Z_{taqif}$ denotes the vector of $K$ explanatory variables describing the
characteristics for this flat (surface in $m^{2}$, count data as number of
rooms, bedrooms, bathrooms, garage plots, and dummy variables as balcony,
whether it is located on a street, boulevard, avenue, or place, period of
construction ($<$1850, 1850-1913, ...,1981-2003), etc). This unbalanced
panel is made up of $N=20$ top-level \textit{arrondissements}, each
containing $Q_{ta}$ second-level \textit{quartiers}. The second-level
\textit{quartiers} in turn contain $M_{taq}$ third-level \textit{\^{\i}lots}%
, which contain the innermost $F_{taqi}$ observations on flats. The number
of observations in the higher level groups are $F_{taq}=%
\sum_{i=1}^{M_{taq}}{}F_{taqi}$ and $F_{ta}=\sum_{q=1}^{Q_{ta}}F_{taq}$.\
The total number of observations is $H=\sum_{t=1}^{T}\sum_{a=1}^{N}F_{ta}$.
The number of top-level groups is $NT,$ the number of second-level groups is
$L=\sum_{t=1}^{T}\sum_{a=1}^{N}Q_{ta}$ and the number of bottom-level groups
is $G=\sum_{t=1}^{T}\sum_{a=1}^{N}\sum_{q=1}^{Q_{ta}}M_{taq}$. So, we have a
\index{five-dimensional}five-dimensional \index{pseudo-panel}pseudo-panel data structure. The \index{spatial!lag}spatial lag coefficient $%
\lambda _{t}$ may be time varying or constant over time and the \index{spatial!lag}spatial lag
variable $\tilde{p}_{taqif}$ is defined as
\begin{equation*}
\tilde{p}_{taqif}=\sum_{a=1}^{N}\sum_{q=1}^{Q_{ta}}\sum_{i=1}^{M_{taq}}%
\sum_{p=1}^{F_{taqi}}w_{taqip}p_{taqip}\,,  \label{CH12X_eq_4a1}
\end{equation*}%
where $w_{taqip}$ denotes the elements of the \index{spatial!time-varying weight}spatial weights
matrices $W_{t}$ which vary with $t$. Elements on the diagonal of $W_{t}$
are set to zero while the off-diagonal elements define the connexion
(\index{contiguity}contiguity or distances) between \index{dwelling}dwellings. There are at least two reasons
why positive spatial correlation may exist. First, \index{dwelling}dwellings in a
\index{neighborhood}neighborhood tend to have similar structural characteristics and second,
\index{dwelling}dwellings in a \index{neighborhood}neighborhood share the same location \index{amenities}amenities (see \citeA
{Basu1988}). However, many of the price determining factors shared by
\index{neighborhood}neighborhoods are difficult to explain explicitly but these
\textquotedblleft omitted" factors are contained in the \index{neighborhood}neighborhood prices.
For each year, \citeA{Baltagi2015}, using the \textquotedblleft Delaunay
triangle algorithm", define first-order \index{contiguity}contiguity matrices $W_{t}$ for the
nearest \index{neighbor}neighbors (\textit{i.e.}, from $10$ to $140$ nearest sold flats).
Accordingly to the \index{nested!structure}nested structure, the disturbance term is given by
\begin{equation*}
\varepsilon _{taqif}=\delta _{ta}+\mu _{taq}+\nu _{taqi}+u_{taqif}\,,
\label{CH12X_eq_4a}
\end{equation*}%
where $\delta _{ta}$ is the \textit{arrondissement} effect, $\mu _{taq}$ is
the \textit{quartier} effect naturally \index{nested!structure}nested in the respective \textit{%
arrondissement} and $\nu _{taqi}$ is the \textit{\^{\i}lot} effect naturally
\index{nested!structure}nested in the respective \textit{quartier}. These could be fixed or random.
The remainder disturbance term for the particular flat is random $%
u_{taqif}\sim iiN(0,\sigma _{u}^{2})$. For the random specification, we
assume that $\delta _{ta}\sim iiN(0,\sigma _{\delta }^{2})$, $\mu _{taq}\sim
iiN(0,\sigma _{\mu }^{2})$ and $\nu _{taqi}\sim iiN(0,\sigma _{\nu }^{2})$.

Following \citeA{Antweiler2001}, \citeA{Baltagi2015} use block-diagonal
matrices of size $\left( H\times H\right) $ corresponding in structure to
the groups or subgroups they represent. They can be constructed explicitly
by using \textquotedblleft group membership\textquotedblright\ matrices
consisting of ones and zeros that uniquely assign each of the $H$
observations to one of the $G$ (or $L$ or $NT$) groups. Let $R_{\nu }$ be
such an $\left( H\times G\right) $ matrix corresponding to the innermost
group level. Then the block-diagonal $\left( H\times H\right) $ matrix $%
J_{\nu }$ can be expressed as the outer product of its membership matrices: $%
J_{\nu }=R_{\nu }R_{\nu }^{\prime }$. The inner product $R_{\nu }^{\prime
}R_{\nu }$ produces a diagonal matrix $\widetilde{L}_{\nu }$ of size $\left(
G\times G\right) $ which contains the number of observations of each group.
Similarly, let $R_{\mu }$ be such an $\left( H\times L\right) $ matrix
corresponding to the second-level groups. Then the block-diagonal $\left(
H\times H\right) $ matrix $J_{\mu }$ can be expressed as the outer product
of its membership matrices: $J_{\mu }=R_{\mu }R_{\mu }^{\prime }$. Last, let
$R_{\delta }$ be such an $\left( H\times NT\right) $ matrix corresponding to
the top-level groups. Then the block-diagonal $\left( H\times H\right) $
matrix $J_{\delta }$ can be expressed as the outer product of its membership
matrices: $J_{\delta }=R_{\delta }R_{\delta }^{\prime }$.\newline
If we pool the observations, the log-likelihood is given by\newline
\begin{equation*}
\ln l=-\frac{1}{2}H\ln \left( 2\pi \right) -\frac{1}{2}\ln |\Omega |+\ln |A|-%
\frac{1}{2}\varepsilon^{\prime }\Omega ^{-1}\varepsilon \,,
\label{CH12X_eq_4b}
\end{equation*}%
where
\begin{equation*}
\varepsilon=Ay-X\beta \text{ , }A=I_{H}-\lambda W \,,  \label{CH12X_eq_4c}
\end{equation*}%
with $W=diag \left( W_{t} \right)$ and $\lambda = diag \left(\lambda _{t}
\right)$ where $W$ is the block-diagonal \index{spatial!weight}spatial weight matrix of size $%
\left( H\times H\right)$. $W_{t}$ is the \index{spatial!time-varying weight}spatial weight matrix\footnote{\citeA
{Baltagi2015} use a block-diagonal weight matrix $W$ of $(156,896\times
156,896)$ whose smallest sub-block is a weight matrix $W_{t}$ of $%
(6,643\times 6,643)$ for the year $1992$ and whose largest sub-block is a
weight matrix $W_{t}$ of $(17,098\times 17,098)$ for the year $1999$.} of
size $\left( F_{ta}\times F_{ta}\right) $ changing at each time period $t$. $%
\lambda $ is the \index{spatial!lag}spatial lag matrix of size $\left( T\times T\right) $ whose
elements $\lambda _{t}$ change at each time
period $t.$ $I_{H}$ is an identity matrix of size $\left( H\times H\right) $.%
\newline
The variance-covariance matrix of the disturbance is defined as follows:
\begin{equation*}
\Omega =E\left[ \varepsilon \varepsilon^{\prime }\right] =\sigma _{u}^{2}%
\left[ I_{H}+\rho _{\nu }J_{\nu }+\rho _{\mu }J_{\mu }+\rho _{\delta
}J_{\delta }\right] \,,  \label{CH12X_eq_4d}
\end{equation*}%
with
\begin{equation*}
\rho _{\delta }=\frac{\sigma _{\delta }^{2}}{\sigma _{u}^{2}}\text{ , }\rho
_{\mu }=\frac{\sigma _{\mu }^{2}}{\sigma _{u}^{2}}\text{ , }\rho _{\nu }=%
\frac{\sigma _{\nu }^{2}}{\sigma _{u}^{2}} \,.  \label{CH12X_eq_4e}
\end{equation*}%
Extending the derivations of \citeA{Antweiler2001} to the case of the spatial
lag model (\ref{CH12X_eq_4}), \citeA{Baltagi2015} get:
\begin{equation}
\ln l=-\frac{1}{2}\left[ H\ln \left( 2\pi \sigma _{u}^{2}\right)
+\sum\limits_{t=1}^{T}\left\{ \ln |I_{t}-\lambda
_{t}W_{t}|+\sum\limits_{a=1}^{N}\left\{ \ln \theta _{ta}+C_{ta}-\frac{\rho
_{\delta }}{\theta _{ta}}\frac{U_{ta}^{2}}{\sigma _{u}^{2}}\right\} \right\} %
\right] \,,  \label{CH12X_eq_4f}
\end{equation}%
\begin{eqnarray*}
\text{with }C_{ta} &=&\sum\limits_{q=1}^{Q_{ta}}\left\{ \ln \theta
_{taq}+C_{taq}-\frac{\rho _{\mu }}{\theta _{taq}}\frac{U_{taq}^{2}}{\sigma
_{u}^{2}}\right\} \,,  \label{equa.7a} \\
\text{and }C_{taq} &=&\sum\limits_{i=1}^{M_{taq}}\left\{ \ln \theta _{taqi}+%
\frac{V_{taqi}}{\sigma _{u }^{2}}-\frac{\rho _{\nu }}{\theta _{taqi}}\frac{%
U_{taqi}^{2}}{\sigma _{u}^{2}}\right\} \,,  \label{CH12X_eq_4g}
\end{eqnarray*}%
where $I_{t}$ is an identity matrix of size $\left( F_{ta}\times
F_{ta}\right) $ and where
\begin{equation*}
\left\{
\begin{array}{ccccccc}
\theta _{taqi} & = & 1+\rho _{\nu }F_{taqi} & \hspace{0.5cm} & V_{taqi} & =
& \sum\limits_{f=1}^{F_{taqi}} \varepsilon_{taqif}^{2} \,, \\
\theta _{taq} & = & 1+\rho _{\mu }\phi _{taq}\text{ } & \text{with }\phi
_{taq}=\left( \sum\limits_{i=1}^{M_{taq}}\frac{F_{taqi}}{\theta _{taqi}}%
\right) \hspace{0.5cm} & U_{taqi} & = & \sum\limits_{f=1}^{F_{taqi}}
\varepsilon_{taqif} \,, \\
\theta _{ta} & = & 1+\rho _{\delta }\phi _{ta}\text{ } & \text{with }\phi
_{ta}=\left( \sum\limits_{q=1}^{Q_{ta}}\frac{\phi _{taq}}{\theta _{taq}}%
\right) \hspace{0.5cm} & U_{taq} & = & \sum\limits_{i=1}^{M_{taq}}\frac{%
U_{taqi}}{\theta _{taqi}} \,, \\
&  &  &  & U_{ta} & = & \sum\limits_{q=1}^{Q_{ta}}\frac{U_{taq}}{\theta
_{taq}} \,,%
\end{array}%
\right.  \label{CH12X_eq_4h}
\end{equation*}%
where $\varepsilon_{taqif}=y_{taqif}-\lambda
_{t}\sum_{a=1}^{N}\sum_{q=1}^{Q_{ta}}\sum_{i=1}^{M_{taq}}%
\sum_{p=1}^{F_{taqi}}w_{taqip}y_{taqip}-X_{taqif}\beta $.\newline
A gradient of this log-likelihood function (\ref{CH12X_eq_4f}) is obtained
analytically, but it can also be obtained through numeric approximation. In
carrying out this maximization, it is necessary to constrain the
optimization such that $|\lambda _{t}|<1,$ the variance $\sigma _{u}^{2}$
remains positive, and that the variance ratios $\rho _{\delta }$, $\rho
_{\mu }$ and $\rho _{\nu }$ remain non-negative.

\citeA{Baltagi2015} report several ML estimation results. One for the \index{random!effects}random effects (RE) model ignoring the \index{nested!effects}nested effects, one for
the \index{nested!random effects model}nested
RE model ignoring the \index{spatial!lag}spatial lag effects and one for the \index{spatial!nested random effects model}spatial nested RE
model.\footnote{%
For the estimation of \index{nested!random effects model}nested error component model with unbalanced panel
data using simple analysis of variance (ANOVA), \index{maximum likelihood}maximum likelihood (MLE) and
minimum norm quadratic unbiased estimators (MINQUE)-type estimators of the
variance components, see \citeA{Baltagi2001}. For Lagrange multiplier testing
of \index{nested!random effects model}nested error component model with unbalanced panel data, see \citeA
{Baltagi2002}.} \citeA{Baltagi2015} found significant \index{spatial!lag}spatial lag effects as
well as significant \index{nested!random effects model}nested random error effects. They emphasize the
importance of \index{nested!effects}nested effects in the \index{housing!model example!Paris}Paris housing data as well as the
\index{spatial!lag}spatial lag effects. In fact, they show that the impact of the adjacent
\index{neighborhood}neighborhoods becomes relatively small when one takes care of the \index{nested!random effects model}nested
\index{random!effects}random effects. In addition, due to the unbalanced \index{pseudo-panel}pseudo-panel aspect of
these transactions, they show that one should allow the \index{spatial!weight}spatial weight
matrix as well as the \index{spatial!lag}spatial lag coefficients to vary over time, and that
the likelihood ratio tests confirm that they fit the \index{housing!model example!Paris}Paris housing data
better.

Following \citeA{LeSage2009}, \citeA{Baltagi2015} compute the \index{marginal effects}marginal effects -- which are decomposed into direct, indirect and total
\index{marginal effects}marginal effects -- and show that the marginal \index{spillover effects}spillover effects due to the
\index{neighbor}neighbors are negligible relative to the \index{direct effects}direct effects. Moreover,
the empirical results show that the \index{marginal effects}marginal effect for a specific housing
characteristic is lower on average once the \index{nested!effects}nested effects are taken into
account.

\citeA{Baltagi2014} \ estimate a \index{nested!random effects model}nested random effects \index{spatial!autoregressive}spatial autoregressive
panel data model to explain annual \index{housing!price}house price variation across 353 local
authority districts in\ England over the period 2000-2007. The \index{nested!structure}nested error
components represent the impact of local (district) effects embedded within
wider (county) effects. \citeA{Baltagi2014} propose new estimators based on
the \index{instrumental variable}instrumental variable approaches of \citeA{Kelejian1998} and \citeA
{Lee2003} for the cross-sectional \index{spatial!autoregressive}spatial autoregressive model. The
estimation methods allow for the \index{endogeneity}endogeneity of the \index{spatial!lag}spatial lag variable
producing the simultaneous spatial spillover of prices across districts
together with the \index{nested!random effects model}nested random effects in a panel data setting. Monte Carlo
results show that these estimators perform well relative to alternative
approaches and produce estimates based on real data that are consistent with
the theoretical \index{housing!price}house price model underpinning the reduced form. The
empirical results show that there is a significant \index{spatial!lag}spatial lag term
indicating positive correlation between prices locally and prices in
`nearby' districts and that income within commuting distance has a positive
effect, while the stock of housing has a negative effect, on \index{housing!price}housing price.
Also that the \index{nested!random effects model}nested error components attributable to district and county
effects, like the \index{spatial!lag}spatial lag term are necessary elements in modeling
\index{housing!model example!UK}UK \index{housing!price}house prices.\footnote{\citeA{BaltagiPirotte2014} derive the Best
Linear Unbiased Predictor (BLUP) for a spatial \index{nested!random effects model}nested error components panel
data model. This predictor is useful for panel data applications that
exhibit \index{spatial!dependence}spatial dependence and a \index{nested!structure}nested hierarchical structure. The
predictor allows for unbalancedness in the number of observations in the
\index{nested!structure}nested groups. It could be interesting for forecasting average housing
prices located in a county \index{nested!structure}nested in a state.} \newline

From \index{hedonic!price}hedonic price functions, we can derive \index{price index!temporal}temporal and/or \index{price index!spatial}spatial price
indexes. This has been done, for instance, by \citeA{Syed2008} for the \index{housing!model example!Sydney}Sydney
region. Their data concern 15 regions in \index{housing!model example!Sydney}Sydney on a quarterly basis from
2001 to 2006 from a data set consisting of 418,877 house sales. As $60\%$ of
sales observations are missing for one or more of the core characteristics, they
first use multiple-imputation techniques to fill in the gaps in the data
set, prior to estimating the hedonic model. In a second stage, they specify
and estimate a non-\index{nested!structure}nested \index{three-dimensional}three-dimensional \index{hedonic!price}hedonic price function. They
pool across all the regions and periods in the sample and estimate the
region-time specific \index{fixed effects}fixed effects and \index{shadow price}shadow prices of housing
characteristics. This method was first proposed by \citeA{Aizcorbe2004}, who
refer to it as the time-interaction-country product dummy method.
\begin{eqnarray*}
p_{jth} &=&\alpha +\sum_{\tau =2}^{T}\beta _{\tau }q_{\tau h}+\sum_{\kappa
=2}^{J}\gamma _{\kappa }r_{\kappa h}+\sum_{\tau =2}^{T}\sum_{\kappa
=2}^{J}\delta _{\tau \kappa }b_{\tau \kappa h}  \notag \\
&&+\sum_{m=2}^{M_{\kappa }}\eta _{\kappa m}d_{\kappa mh}+Z_{jth}\theta
+\varepsilon _{jth}\,,  \label{CH12X_eq_5} \\
&&\text{ for }j=1,...,J\text{, }t=1,...,T\text{ and }h=1,...,H_{jt}\,,
\notag
\end{eqnarray*}%
where $p$ is the log of the price of a \index{dwelling}dwelling $h$ belonging to
region-period $jt$, $q_{\tau h}$ (resp. $r_{\kappa h}$) are dummy variables
such that $q_{\tau h}=1$ (resp. $r_{\kappa h}=1$) if the observation $h$ is
from period $t$ (resp. from region $j$) and zero otherwise. The dummy
variables $b_{\tau \kappa h}$ denote interactions between periods and
regions taking the value of $1$ if the observation $h$ is from region-period
$jt$ and zero otherwise. The postcode dummies are denoted by $d_{\kappa mh}$%
, where $d_{\kappa mh}=1$ for observation $h$'s postcode and zero otherwise.
$Z$ is a set of \index{quality}quality characteristics including the \index{dwelling}dwelling type, the
number of bedrooms, bathrooms, lot size, etc. Spatial correlation between
observations is defined by a \index{spatial!autoregressive}spatial autoregressive process on the error
term: $\varepsilon _{jth}=\lambda W\varepsilon _{jth}+u_{jth}$ where $%
u_{jth}\sim N(0,\omega _{jth}\sigma ^{2})$. The \index{spatial!weight}spatial weight matrix $W$ is
a \index{contiguity}contiguity matrix and the variance of $u_{jth}$ is subscripted with $jt$
allowing for heteroskedasticity. The coefficients $\delta _{jt}$ measure the
region-period specific \index{fixed effects}fixed effects for the logarithms of the price level
after controlling for the effects of the \index{attributes}attributes of the \index{dwelling}dwellings. The
model is estimated using the \index{maximum likelihood}maximum likelihood method. The advantage of
this region-time-dummy model is that the \index{price index!temporal}temporal and \index{price index!regional}regional price indexes are derived directly from the estimated coefficients $\hat{\beta}%
_{t}$, $\hat{\gamma}_{j}$, $\hat{\delta}_{jt}$, $\hat{\eta}_{jm}$ and $\hat{%
\theta}$. Let $P_{j,t,s}$ the price index for region $j$ in year $t$ and
quarter $s$. Then, the relative prices are given by
\begin{eqnarray*}
\frac{P_{j,t,s}}{P_{j,t,1}} &=&\exp \left( \hat{\beta}_{t,s}+\hat{\delta}%
_{j,t,s}\right) \text{ for }s=2,3,4\,,  \notag \\
\text{ and }\frac{P_{j,t+1,1}}{P_{j,t,1}} &=&\exp \left( \hat{\beta}_{t+1,1}+%
\hat{\delta}_{j,t+1,1}\right) \,.  \label{CH12X_eq_6}
\end{eqnarray*}%
So, it is possible to construct a \index{price index!temporal}temporal price index for each
region $j$ over the entire time period of the dataset. Results are
normalized such that the price index for the initial region (Inner \index{housing!model example!Sydney}Sydney)
is equal to $1$ for the first quarter of 2001. One can also construct a
\index{price index!spatial}spatial price index for each quarter $s$ of a specific year $t$ for the
entire set of regions. For a given quarter ($t,s$), spatial price indexes
can be constructed from the estimated coefficients $\hat{\gamma}_{j}$, $\hat{%
\delta}_{jt}$, $\hat{\eta}_{jm}$ and $\hat{\theta}$. The starting point is a
\index{housing!model example!Paris}comparison between a postcode $m$ in region $l$ and a postcode $n$ in region
$j$ for a particular \index{dwelling}dwelling $h$ with \index{amenities}amenities vector $Z_{ch}$. This
\index{price index!spatial}spatial price index is defined as:
\begin{eqnarray*}
P_{lmts,jnts}(Z_{ch}) &=&\exp \left[ \left( \hat{\gamma}_{j}-\hat{\gamma}%
_{l}\right) +\left( \hat{\delta}_{jt}-\hat{\delta}_{lt}\right) +\left( \hat{%
\eta}_{jn}-\hat{\eta}_{lm}\right) \right]   \notag \\
&&\times \left[ \prod_{c=1}^{C}\exp \left[ Z_{ch}\left( \hat{\theta}_{jc}-%
\hat{\theta}_{lc}\right) \right] \right] \,,  \label{CH12X_eq_7}
\end{eqnarray*}%
and the spatial index can be generalized to take into account of all
\index{dwelling}dwellings sold in postcodes $lm$:
\begin{eqnarray*}
P_{lmts,jnts} &=&\exp \left[ \left( \hat{\gamma}_{j}-\hat{\gamma}_{l}\right)
+\left( \hat{\delta}_{jt}-\hat{\delta}_{lt}\right) +\left( \hat{\eta}_{jn}-%
\hat{\eta}_{lm}\right) \right]   \notag \\
&&\times \left[ \prod_{h=1}^{H_{lmts}}\prod_{c=1}^{C}\exp \left[
Z_{ch}\left( \hat{\theta}_{jc}-\hat{\theta}_{lc}\right) \right] \right]
^{1/H_{lmts}}\,.  \label{CH12X_eq_8}
\end{eqnarray*}%
This is close to a \index{price index!Laspeyres}.

Combining the \index{price index!temporal}temporal and \index{price index!spatial}spatial indexes allows a price \index{housing!model example!Paris}comparison of
\index{dwelling}dwellings between different location-year-quarter triplets. \citeA{Syed2008}
found that their hedonic \index{housing!price}house price indexes rose significantly from 2001 to
2003, after which they fell slightly. This finding is consistent with the
Australian Bureau of Statistics (ABS) index. But their indexes, however, are
less volatile than their ABS counterpart, rising noticeably less in the boom
and falling less thereafter. In the spatial dimension, they found large and
systematic differences in the price of housing across regions of \index{housing!model example!Sydney}Sydney. The
regional dispersion narrowed during the boom period but appears to have
increased again since then.\vspace{0.5cm}

Several authors have shown that values of complex assets are difficult to
accurately quantify and \index{information asymmetry}information asymmetry affects asset prices
through various channels (see \citeA{Agarwal2010}, \citeA{Baker2007}, \citeA
{Carlin2013}, \citeA{Kelly2012} to mention a few). The \index{subprime crisis}subprime crisis %
(poor household \index{mortages}mortgage decisions and subsequent foreclosure), the
housing market collapse in the US, followed by the \index{financial crisis}financial crisis have
revealed that uninformed buyers overpay. The house buying mechanism is a
field in which households' ability (or inability) to use market information
may have strong effects on housing decisions. It could be done through the
choice of \index{mortages}mortgage product and through the purchase transaction (see \citeA
{Carlin2013} and \citeA{Turnbull2015}). House purchases may involve
residential \index{mortages}mortgages and associated complex financial instruments which
have been assigned as a major cause of waves of \index{foreclosures}foreclosures during
and after the 2007-2008 \index{financial crisis}financial crisis. \citeA{Turnbull2015} show
that buyers who are uninformed in the housing market pay more for houses
than buyers who are informed. They use \index{pseudo-panel}pseudo-panels of \index{repeated sales}repeated sales based on \index{neighborhood}neighborhood
census block-level. This data is for 426,021
parcels located in \index{housing!model example!US Florida}Orange County, Florida, over the period 2000-2012. The
authors split fair market value and uninformed buyer effects by first
identifying for each of the market sales in the period 2000-2006 which of
the units foreclosed in 2007-2012. The future foreclosure dummy $FF$ equals $%
1$ if a market transaction completed in 2000-2006 is followed by a
foreclosure in 2007-2012 and equals zero otherwise. \citeA{Turnbull2015}
estimate an \index{hedonic!price}hedonic price function of the log of market price in first
differences on the \index{neighborhood}neighborhood block-level $j$:
\begin{eqnarray*}
p_{itj}-p_{lsj} &=&\left( Z_{itj}-Z_{lsj}\right) \beta _{Z}+\left(
FF_{ij}-FF_{lj}\right) \beta _{FF}+\varepsilon _{itj}-\varepsilon _{lsj}\,,
\label{CH12X_eq_9} \\
&&\text{ for }t,s=1,...,T\text{, }i,l=1,...,N\text{, }j=1,...,J\text{ for
all }i\neq j\text{ and }t\neq s\,,  \notag
\end{eqnarray*}%
where $p_{itj}$ is the log of the price of property $i$ sold at time $t$
located in area $j$. $Z$ is the vector of relevant house characteristics and
\index{amenities}amenities and $FF$ is the penalty associated with being \index{foreclosures}foreclosed \textit{%
ex post} (over 2007-2012). The model of first differences at the
\index{neighborhood}neighborhood block-level basically treats sales within the \index{neighborhood}neighborhood
block as \index{repeated sales}repeat sales while accounting for observed structural differences.
This is a model on \index{pseudo-panel}pseudo-panels of repeated observations \textquotedblleft
\textit{\`{a} la} \citeA{Deaton1985}". This model also allows for \index{clustered errors}clustered errors at the \index{neighborhood}neighborhood block-level $j$. Results show that
buyers who are later \index{foreclosures}foreclosed paid a $2.7\%$ (resp. a $4.6\%$) premium for
properties bought between 2000 and 2006 (resp. between 2005 and 2006).
Estimation on different sub-periods also reveal strong correlation between
home buyers' \index{housing!price}house prices and future \index{foreclosures}foreclosures. To check whether effects
vary across housing market segments, \citeA{Turnbull2015} estimate \index{quantile regression}quantile regression models. Results show that the effect for the penalty associated
with being foreclosed is larger for the lower end of the housing market.
Buyers in 2005-2006 who ended up foreclosed paid up to $3.5\%$ above fair
market value in the lower end of the housing market while foreclosed owners
paid a little over $1\%$ percent more in the higher end of the housing
market.\\


Evaluating non-marginal changes in amenities in the hedonic model requires estimating the underlying hedonic demand or \index{marginal willingness-to-pay function}marginal willingness-to-pay (MWTP) function. In general, however, it is necessary to make assumptions that restrict preference heterogeneity. As stated earlier, the usual approach, based on \citeA{Rosen1974}'s seminal paper, suffers from a number of econometric problems, especially when the hedonic price function is non-linear in the amenity of interest or when buyers simultaneously choose both the hedonic price and the quantity of that amenity that they will consume. The introduction of a nonlinear budget constraint creates a serious endogeneity problem when using statistical inference to recover the parameters describing the consumer's preferences. And the instrumental-variable approaches generally used have relied on questionable exclusion restrictions. These problems have led most researchers to abandon the estimation of the MWTP function altogether and to use only local MWTP measurements. 
But, \citeA{bajari2005demand} have demonstrated that this endogeneity problem may be avoided by replacing the statistical inference used in the second stage of \citeA{Rosen1974}'s method with a ``preference inversion'' procedure that inverts the first-order conditions of utility maximization to recover demand at the individual level. However, strict functional-form assumptions, necessary for their inversion procedure, limit their approach. In contrast, \citeA{bishop2018using} have shown how to recover the unconditional distribution of linear MWTP functions with a simple and transparent data-driven estimation approach. Constructing a rich 3D panel dataset of property transactions and of buyers covering the six core counties of the \index{housing!model example!US San Francisco}San Francisco Bay Area, they recover the full distribution of demand functions for clean air. First, they find that estimating the full demand function, rather than simply recovering a local estimate of marginal willingness to pay, is important. Second, they find  that the mean willingness to pay to avoid a $33\%$ increase in ozone is $\$1,021$ and that the median willingness to pay to avoid the increase is $\$770$. There is considerable heterogeneity as evidenced by the interquartile range of $\$1,106$. Their data-driven estimates of the welfare effects associated with a nonmarginal change in air quality differ substantially from those recovered using the existing approaches to welfare estimation.



\section{Multi-dimensional Models of Residential Mobility and Location
Choice: Some Examples}

\label{CH12X_section4}

\index{residential!mobility}Residential mobility and \index{location choice}location choice are significant
household decisions and have been widely researched in various fields
including economics, sociology, geography, regional science, urban planning,
housing policy, transportation, etc. Decisions of \index{residential!mobility}residential mobility and
\index{location choice}location choice are closely related to the household housing process with a
large range of factors that contribute to each choice. Due to the vastness
of the literature on such topics, we will focus on few examples of
\index{residential!mobility}residential mobility and \index{location choice}location choice. Readers could read profitably the
survey of \citeA{Dieleman2001} on \index{residential!mobility}residential mobility. Since the seminal
works of \citeA{Rossi1955} and \citeA{Alonso1964}, a huge number of research
on residential \index{location choice}location choice has been published. \textit{\textquotedblleft
Reasons for moving are divided into those which pertain to the decision to
move out of the former home - \textquotedblleft pushes" - and those reasons
pertaining to the choice among places to move to - \textquotedblleft pulls" }
(\citeA{Rossi1955}, p. 8). For instance, push factors may include negative
externalities like noise, pollution or crime, changes in housing
affordability, dissatisfaction with the current \index{dwelling}dwelling, changes in
household structure, etc. Pull factors often include better access to good
\index{quality}quality public services (schools and health care facilities), employment,
leisure and recreational opportunities, etc. (see \citeA{Lee2010} and \citeA
{Hoang2000} for a review). Our purpose is not to review the main factors of
\index{residential!mobility}residential mobility and relocation but to summarize few multi-dimensional
studies of \index{residential!mobility}residential mobility and relocation.

One interesting study has been done by \citeA{Davies1985} in a
multi-dimensional framework. They propose a model that conceptualizes
\index{residential!mobility}residential mobility as a sequence of choices between staying and moving.
Household $i$ will move in time period $t$ if and only if \index{random!utility}random utility
derived from the most-favored alternative \index{dwelling}dwelling available $u_{itb}$ is
larger than the \index{random!utility}random utility derived from the current \index{dwelling}dwelling $u_{ita}$:
\begin{eqnarray*}
u_{ita} &=&V\left( y_{it},Z_{ta}\right) +\varepsilon
_{ita}=V_{ita}+\varepsilon _{ita}\text{ with }\varepsilon _{ita}=\mu
_{ia}+g(d_{it})+\nu _{ita} \,,  \label{CH12X_eq_10} \\
u_{itb} &=&V\left( y_{it},Z_{tb}\right) +\varepsilon
_{itb}=V_{itb}+\varepsilon _{itb}\text{ with }\varepsilon _{itb}=\mu
_{ib}+h(t)+\nu _{itb} \,,  \notag
\end{eqnarray*}%
where $y_{it}$ is a vector of observed characteristics of household $i$ at
time $t$, $Z_{ta}$ (resp. $Z_{tb}$) is a vector of observed characteristics
of the current \index{dwelling}dwelling (resp. the most-favored alternative \index{dwelling}dwelling
available). $V_{ita}$ and $V_{itb}$ are the systematic utilities while $%
\varepsilon _{ita} $ and $\varepsilon _{itb}$ are the random components of
utilities. These random components are likely to be correlated over time for
each household. $\varepsilon _{ita}$ is the sum of the unexplained household
\index{heterogeneity}heterogeneity $\mu _{ia}$, a function $g(d_{it})$ of the \index{duration!of stay}duration of stay
for household $i$ at time $t$ and a remainder term $\nu _{ita}$,
independently distributed over both households and time. For the other
random component $\varepsilon _{itb}$, the unexplained household
\index{heterogeneity}heterogeneity $\mu _{ib}$ also applies. Moreover, a time trend $h(t)$
represent fluctuations in market conditions. \citeA{Davies1985} used a
quadratic specification for the \index{duration!of stay}duration of stay $g(d_{it})=\beta
_{1}d_{it}+\beta _{2}d_{it}^{2}$, and a cubic specification for the housing
market function $h(t)=\beta _{3}t+\beta _{4}t^{2}+\beta _{5}t^{3}$.

The likelihood $L(z_{it})$ of the observed sequence of outcomes is the
product of the probabilities of the observed choice for each time period:
\begin{eqnarray*}
L(z_{it}) &=&\prod_{t=1}^{T}\left\{ Pr\left[ u_{itb}>u_{ita}\right] \right\}
^{z_{it}}\left\{ 1-Pr\left[ u_{itb}>u_{ita}\right] \right\} ^{1-z_{it}}\,,
\label{CH12X_eq_11} \\
\text{with } &&Pr\left[ u_{itb}>u_{ita}\right] =\int_{-V_{it}-\mu
_{i}+g(d_{it})-h(t)}^{\infty }\phi \left( \nu _{itb}-\nu _{ita}\right)
d\left( \nu _{itb}-\nu _{ita}\right) \,,  \notag
\end{eqnarray*}%
where $z_{it}=1$ if household $i$ moves in time period $t$ and zero
elsewhere, $V_{it}=V_{itb}-V_{ita}$, $\mu _{i}=\mu _{ib}-\mu _{ia}$ and $%
\phi \left( .\right) $ is the probability density of the difference between
the two random components. Assuming that they follow \index{Weibull distribution}Weibull distributions leads to the following likelihood with a household-specific
error term $\mu _{i}$:
\begin{equation*}
L(z_{it})=\prod_{t=1}^{T}\frac{\exp \left[ -V_{it}-\mu _{i}+g(d_{it})-h(t)%
\right] ^{z_{it}}}{1+\exp \left[ -V_{it}-\mu _{i}+g(d_{it})-h(t)\right] }\,.
\label{CH12X_eq_11a}
\end{equation*}%
Three problems arise with this likelihood: the integration over the error
term distribution is almost analytically intractable; the initial
observation complicates the handling of endogenous variables such as
\index{duration!of stay}duration of stay $d_{it}$ and numerical methods are required for parameter
estimation. To overcome these problems, \citeA{Davies1985} derived an
approximation of the likelihood using the \index{generalized!Beta-logistic approach}generalized Beta-logistic approach
developed by \citeA{Davies1984}.

The panel data is for 887 households participating in the \index{Michigan Panel Study of Income Dynamics}Michigan Panel Study of Income Dynamics over the period 1968-1977.
The dependent
variable was a residential move within the county or an intercounty move
with no change in the head-of-household's job. Among the main explanatory
variables were the \index{duration!of stay}duration of stay, a room adequacy index (actual rooms /
required rooms), an income adequacy index (actual income / needs), the age
of the head of household, and the education level. First, they show that the
room adequacy index has a U-shaped relationship with \index{residential!mobility}residential mobility.
Renters have the shortest initial duration status while owners have the
longest. But, there is no evidence of a similar U-shaped relationship
anticipated for the income adequacy index. Second, they show that changing
financial circumstances does not seem to play any role in the life cycle
variation in \index{residential!mobility}residential mobility in the United States. Moreover, they are
not able to show any effect of income surplus on \index{residential!mobility}residential mobility. These
are unexpected results. \citeA{Davies1985} argue that these results may be
due to the housing market being highly segmented, not just between renting
and owner-occupation, but between different types of property and their
location. It could be interesting to redo this study on more recent data. It
will probably give different conclusions for the last decade which has known
troubled financial periods. \citeA{Davies1985} found a strong negative
relationship between age of the head of household and \index{residential!mobility}residential mobility.
This strong negative relationship is present even when changing space
requirements and financial pressures are accounted for. Age of the head of
household is the dominant life cycle and acts as a proxy variable for
changing needs and financial circumstances through the life cycle.

Several studies have acknowledged the interdependencies between residential
choice, job choice, and transportation mode choice. However, most previous studies focused on two dimensions only. \citeA{guo2020co} jointly consider a multidimensional model of residential choice, job choice, commuting mode choice and assume individuals and households consider different life domains jointly. Based on a 3D dataset (individuals, choice situations and alternatives) collected in the region of Shenyang (China), \citeA{guo2020co} specify a \index{mixed multinomial logit model}mixed multinomial logit model which allows for unobserved heterogeneity in individual preferences. This model was estimated to capture the effects of different residential, job and commuting attributes on multi-dimensional choice, accounting for the panel nature of the data. \citeA{guo2020co} find that, housing tenure, size, price, distance to the bus stop, and housing location are important housing characteristics explaining the residential mobility choice process. Second, salary, job type, co-worker relationships and job environment are also significant factors in the job mobility choice process. \citeA{guo2020co} show that time-related factors influence commute mode choice and choice of public transportation modes is sensitive to commuting costs while car mode choice is not. The model estimation highlights that people are relatively satisfied with their current situation and do not frequently make changes. And people are less inclined to move house relative to changing job. Finally, \citeA{guo2020co} show that both unobserved heterogeneity and demographic characteristics affect the multiple dimensions of choice.

Explaining the factors which determine housing tenure choices is important.
For instance, \citeA{Fu2015} estimate \index{multilevel multinomial logistic regression}multilevel multinomial logistic regressions for housing types to study home ownership in urban
\index{housing!model example!China}China. They
base their estimation on a sample of $2,585,480$ households from the
2005 \index{National Population Sample Survey of China}National Population Sample Survey of China and available
information for 205 urban areas (prefectures-level data) (see \citeA
{Huang2002} for a similar study in China but in a 2D framework). For one
household $i$ in prefecture $j$, the within-prefecture multinomial logistic
model for the odds of housing type $m$ is given by
\begin{eqnarray*}
\log \left[ \frac{Pr\left( \text{housing type}_{mij}\right) }{Pr\left( \text{%
private rental housing}_{ij}\right) }\right] =\beta
_{mj0}+\sum_{k=1}^{K}\beta _{km}\left( Z_{h,kijm}-\overline{Z}%
_{h,kjm}\right) +\varepsilon _{ij}\,. &&  \notag \\
&&  \label{CH12X_eq_12a}
\end{eqnarray*}%
The $m=1,..,5$ housing types refer to owning self-built housing, owning
commodity housing, owning affordable housing, owning privatized \textit{%
danwei} housing and public rental housing. $Z_{h,kijm}$ is the value of
household-level covariate $k$ associated with household $i$ in prefecture $j$
for the $m$-th housing type. $\overline{Z}_{h,kjm}$ is the sample mean of
covariate $k$ within prefecture $j$. The household-level error term $%
\varepsilon _{ij}$ is assumed to be $i.i.N(0,\sigma ^{2})$. The
between-prefecture model for housing types is:
\begin{equation*}
\beta _{mj0}=\gamma _{00m}+\sum_{s=1}^{S}\gamma _{0sm}Z_{p,sjm}+\eta
_{0jm}\,,  \label{CH12X_eq_12b}
\end{equation*}%
where $Z_{p,sjm}$ is the prefecture-level covariate $s$ in prefecture $j$
for the $m$-th housing type and $\eta _{0jm}$ is the prefecture-level error
term, which is assumed to be $i.i.N(0,\sigma _{m}^{2})$.

Using \index{generalized!linear mixed model}generalized linear mixed model with \index{random!effects}random effects estimation
methods (GLMM), \citeA{Fu2015} show, at the household level, that
redistributors (\textit{e.g.}, cadres) and supporting clerical staff were
more likely to achieve home ownership than manual workers did. Both
non-agricultural status and working in state sectors confer benefits in
obtaining reform-era housing with heavy subsidies or better qualities. When
one takes into account education and earnings, the advantage of
redistributors (\textit{e.g.}, cadres) over manual workers in home ownership
could be explained by work units. At the prefecture-level, they show that
the marketization only reduced the local home ownership of self-built
housing, affordable housing and privatized \textit{danwei} housing but not
that of commodity housing. In contrast, political and market connections
promote all types of home ownership except self-built housing, and have
significant positive association with the odds of renting public housing.

In the literature, numerous studies focus on how \index{neighborhood}neighborhoods change in
terms of income level, housing values, environment \index{amenities}amenities or different
racial preferences, etc. Racial and ethnic composition may have effects on
\index{neighborhood}neighborhood economic change (see for instance \citeA{Sykes2003}). Some
studies have examined how \index{neighborhood}neighborhood minority composition is associated
with change in \index{neighborhood}neighborhood relative economic status. This is, for instance,
the case of the paper of \citeA{Jun2016} in a 3D framework. He uses the
\index{Neighborhood Change Database}Neighborhood Change Database (NCDB), which includes the decennial
census data across the \index{housing!model example!US}USA from 1970 to 2000 at the census tract level. The
multilevel modeling fits the data structure that a \index{neighborhood}neighborhood is \index{nested!structure}nested in
a metropolitan area and allows for answering the research question whether
the effect of \index{neighborhood}neighborhood racial/ethnic composition on \index{neighborhood}neighborhood
economic change is conditioned by metropolitan-level factors. \citeA{Jun2016}
shows that both \index{neighborhood}neighborhood percentage Black and Hispanic are negatively
related to \index{neighborhood}neighborhood economic gain and are conditioned by
metropolitan-level factors. Although this negative effect of \index{neighborhood}neighborhood
minority composition has been consistent over the four ten years panel, --
the 1970s, 1980s, 1990s, and 2000s -- its impact level is lower in the
latest panel compared to the early panel. The negative effect of
\index{neighborhood}neighborhood minority composition has also declined by the interactions with
metropolitan minority composition. In the later panels, metropolitan
minority composition turned out to moderate the negative effect of
\index{neighborhood}neighborhood minority composition.


We now turn to studies focusing on residential mobility, location choice and the impact of the \index{U.S. housing choice voucher program}U.S. housing choice voucher program.\footnote{The Section 8  housing choice voucher program is the federal government's major program for assisting very low-income families, the elderly, and the disabled to afford decent, safe, and sanitary housing in the private market. Since housing assistance is provided on behalf of the family or individual, participants are able to find their own housing, including single-family homes, townhouses and apartments. The participant is free to choose any housing that meets the requirements of the program and is not limited to units located in subsidized housing projects. A family that is issued a housing voucher is responsible for finding a suitable housing unit of the family's choice where the owner agrees to rent under the program.} \citeA{eriksen2015housing} estimate the effect of increasing the \index{U.S. housing choice voucher program!supply of housing vouchers}supply of housing vouchers on rents using a panel of housing units in the American Housing Survey. Their full-sample is a 3D panel dataset of $8,388$ rental housing units located in the $135$ identified MSAs in the public use version of the AHS with a reported rent in 1997 and at least one other year between 1999 and 2003. \citeA{eriksen2015housing} find an elasticity which is close to zero. When splitting the sample at arbitrary values of a relevant variable (ratio of the rent of a rental unit in the base year to the US Department of Housing and Urban Development fair market rent for that metropolitan statistical area in the same year), they find statistically significant negative effects for those units which were initially $80\%$ below the 1997 ratio and statistically significant positive elasticities for those units which were between $80\%$ and $120\%$ of that ratio. They do not find that an increase in vouchers affected the overall price of rental housing but do estimate differences in effects based on an individual unit's rent before the voucher expansion. Their results are consistent with voucher recipients renting more expensive units after receiving the subsidy. They also find that the largest price increases were for units near the maximum allowable voucher rent in cities with an inelastic housing supply. Several research studies have shown that expanding the housing choice voucher program has varying price effects on different parts of the rental housing market. While recipients will certainly benefit from the added vouchers, there is ambiguity in the impacts on quality of life for non-recipients in the broader housing market. Excessive crowding in homes, dubbed overcrowding, is one measure of quality of life because overcrowding can lead to adverse outcomes such as stunted child development, adverse mental health, and an increased risk in spreading infectious diseases such as COVID-19.\\ 
\citeA{kole2022housing} makes use of an exogenous increase in the \index{U.S. housing choice voucher program!supply of housing vouchers}supply of housing vouchers to explore its effects on overcrowding in the general housing market in the U.S. Using the dataset of \citeA{eriksen2015housing}, estimation of a \index{linear probability model}linear probability model shows that adding vouchers reduces the incidence of overcrowding: a $10\%$ increase in the supply of vouchers reduces the likelihood that a housing unit is overcrowded by $0.081$ percentage points. The mechanism behind the effect is shown to correspond with anecdotes about overcrowding: households experiencing misfortune, financial difficulties, or other tenuous circumstances double-up with higher-income households. Hence, a voucher enables the troubled household to move into more suitable living arrangements.

\citeA{henderson2022nonparametric} develop a nonparametric estimator that works for multidimensional fixed effects model, has a closed form solution and can be estimated in a single step. They apply their method on the relationship between the price of rental housing and housing vouchers using the dataset of \citeA{eriksen2015housing}. As previously stated, an heterogeneous result with respect to the ratio of the rent of a rental unit to the US Department of Housing and Urban Development fair market rent has been observed in the literature (\citeA{eriksen2015housing}) using arbitrary splits of the sample. \citeA{henderson2022nonparametric} avoid these arbitrary splits by adopting a semiparametric approach whereby they obtain an elasticity for each rental unit in the sample and confirm existence of both negative and positive impacts of housing vouchers on the price of rental housing. They find that positive elasticities are concentrated in the Western United States and specifically in areas which are more supply inelastic, but overall, negative elasticities are more prominent in this dataset. This suggests that increasing rents for those who do not receive subsidies are likely localized and not predominant in the U.S. \\


Explaining the \index{residential!choice}residential choices and the \index{residential!mobility}residential mobility is not
sufficient. It seems important to jointly model the \index{residential!mobility}residential mobility and
the \index{duration!of stay}duration of stay at a location preceding the relocation. A lot of
research treated the decision to move as a \index{discrete choice}binary choice decision
(move/no-move) and modeled this decision as a function of various factors
(see above). Others have used \index{duration!model}duration models (see \citeA{Deng2003}) to
represent the stay at a location between moves, treating the reason for a
move as an exogenous variable. An interesting study, done in a
multi-dimensional framework by \citeA{Eluru2009}, has extended these previous
studies in three ways. First, the move decision is treated as an endogenous
variable in a multinomial unordered choice modeling framework. Second, the
\index{duration!of stay}duration of stay is modeled as a grouped choice, supposing that households
treat the \index{duration!of stay}duration of stay at a \index{residential!location}residential location in terms of time-period
ranges as opposed to exact continuous durations. Third, they consider
\index{heterogeneity}heterogeneity of exogenous variables using \index{random!coefficients}random coefficients in both the
equation for the move as well as the equation for the \index{duration!of stay}duration of stay
preceding a relocation. In sum, \citeA{Eluru2009} estimated a joint unordered
choice-grouped choice model system with \index{random!coefficients}random coefficients.

Let the households be represented by the index $i=1,...,N$, let the
different move reasons (\textit{e.g.}, personal reasons, employment reasons,
etc.) be represented by the index $m=1,...,M$ and let the duration
categories (\textit{e.g.}, $<2$ years, $2-5$ years, $5-10$ years, etc.) be
represented by the index $j=1,2,...,J$. The specification of \citeA{Eluru2009}
allows the possibility of multiple move records per household defined by the
index $t=1,2,...,T$ as the different moving choice occasions for households $%
i$. The system of equations jointly models the reason for move and the
\index{duration!of stay}duration of stay as follows:
\begin{equation*}
\left\{
\begin{array}{lll}
u_{imt} & = & X_{it}^{\prime }\beta _{im}+\eta _{im}+\varepsilon _{imt} \,, %
\label{CH12X_eq_13} \\
d_{imt} & = & j \text{ if }\psi _{m,j-1}<d_{imt}^{\ast }<\psi _{m,j} \,,%
\end{array}%
\right.
\end{equation*}

\begin{equation*}
\begin{array}{lll}
\text{ with }d_{imt}^{\ast } & = & X_{it}^{\prime }\alpha _{im}\pm \eta
_{im}+\zeta _{imt} \,. \label{CH12X_eq_13a}%
\end{array}%
\end{equation*}

The first equation of the system is associated with the \index{random!utility}random utility $%
u_{imt}$ for a household $i$ corresponding to the reason to move $m$ at
choice occasion $t$. The $(Q\times 1)$ vector $X_{it}$ is the vector of
\index{attributes}attributes associated with household $i$ and its choice environment (\textit{%
e.g.}, sex, age, employment status, family type, transportation mode to
work, etc.) at the $t$-th choice occasion. The $(Q\times 1)$ random
coefficient vector $\beta _{im}=\beta _{m}+\gamma _{im}$ is the sum of a
vector $\beta _{m}$ of mean effects of the elements of $X_{it}$ for move
reason $m$ and a random vector $\gamma _{im}$ with its $q$-th element ($%
q=1,...,Q$) representing unobserved factors specific to household $i$ and
his choice environment. $\eta _{im}$ expresses unobserved individual factors
that simultaneously impact the propensity of moving for a certain reason $m$
and the \index{duration!of stay}duration of stay. $\varepsilon _{imt}$ is an idiosyncratic random
error term assumed to be identically and independently standard \index{Gumbel distribution}Gumbel distributed across individuals, move reasons and choice occasions.

The second equation of the system is associated with $d_{imt}^{*} $, being
the latent (continuous) \index{duration!of stay}duration of stay for household $i$ before moving for
reason $m$ at the $t$-th choice occasion. This latent duration is mapped to
the grouped duration category $d_{imt}$ by the $\psi$ thresholds (with
infinite bounds as in the usual ordered-response modeling framework). $%
d_{imt}$ is observed only if the end of the \index{duration!of stay}duration of stay at a
\index{residential!location}residential location is associated with alternative $m$. The $(Q \times 1)$
\index{random!coefficients}random coefficient vector $\alpha_{im} = \alpha_{m} + \delta_{im}$ is the
sum of the vector $\alpha_{m}$ of mean effects for category $m$, and the
random vector $\delta_{im}$ of unobserved factors specific to household $i$
and his \index{duration!of stay}duration of stay. $\zeta_{imt}$ is an idiosyncratic random error
term, assumed identically and independently logistic distributed across
individuals, reasons for move, and choice occasions, with variance $%
\lambda^2 $. The elements of the random vectors $\gamma$, $\delta$ and $\eta$
are normally distributed: $\gamma_{imq} \sim N(0, \sigma_{\gamma_{mq}}^2)$, $%
\delta_{imq} \sim N(0, \sigma_{\delta_{mq}}^2)$ and $\eta_{im} \sim N(0,
\sigma_{\eta_{m}}^2)$ for $q=1,...,Q$.

Correlation in unobserved individual factors between the reason to move and
the \index{duration!of stay}duration of stay may be positive or negative, it is indicated by the $%
\pm $ sign in front of $\eta _{im}$ in the duration category equation. If a
positive sign seems logical for the propensity of a move for a given reason $%
m$ in the first equation, a negative sign in the second equation suggests
that unobserved individual factors will decrease the \index{duration!of stay}duration of stay
preceding such a potential move. In the estimation, \citeA{Eluru2009}
considered both the positive and negative signs on the $\eta _{im}$ terms in
the second equation of the system. But the negative sign for all $m$
provided statistically superior results. Conditional on $\gamma _{im}$ and $%
\eta _{im}$ for each (and all) $m$, the probability of a household $i$
choosing to move for reason $m$ at the $t$-th choice occasion is given by
\begin{equation*}
P_{imt}=\frac{\exp \left( X_{it}^{\prime }\beta _{im}+\eta _{im}\right) }{%
\sum_{m=1}^{M}\exp \left( X_{it}^{\prime }\beta _{im}+\eta _{im}\right) }\,.
\label{CH12X_eq_13b}
\end{equation*}%
Conditional on $\delta _{im}$ and $\eta _{im}$, the probability of a
household $i$ choosing to stay for a particular duration category $j$
preceding a move for reason $m$ at the $t$-th choice occasion is given by
\begin{equation*}
R_{imtj}=G\left( \frac{\psi _{m,j}-\left\{ X_{it}^{\prime }\alpha _{im}\pm
\eta _{im}\right\} }{\lambda }\right) -G\left( \frac{\psi _{m,j-1}-\left\{
X_{it}^{\prime }\alpha _{im}\pm \eta _{im}\right\} }{\lambda }\right) \,,
\label{CH12X_eq_13c}
\end{equation*}%
where $G(.)$ is the cumulative distribution of the standard logistic
distribution. Let $\Omega $ be a vector that include all the parameters $%
\beta _{m}$, $\alpha _{m}$, $\lambda $, $\sigma _{\gamma _{mq}}$, $\sigma
_{\delta _{mq}}$ and $\sigma _{\eta _{m}}$ for $m=1,...,M$ and $q=1,...,Q$.
Let $c_{i}$ be a vector stacking the coefficients $\gamma _{im}$, $\delta
_{im}$ and $\eta _{im}$ across all $m$ for household $i$. Let $\Sigma $ be
another vector stacking the standard error terms $\sigma _{\gamma _{mq}}$, $%
\sigma _{\delta _{mq}}$ and $\sigma _{\eta _{m}}$ and let $\Omega _{-\Sigma }
$ represent a vector of all parameters except the standard error terms.
Then, the unconditional likelihood function for all the households is given
by
\begin{eqnarray}
L\left( \Omega \right) =\prod_{i=1}^{N}L_{i}\left( \Omega \right)
&=&\prod_{i=1}^{N}\int_{c_{i}}\left\{ L_{i}\left( \Omega _{-\Sigma }|
c_{i}\right) \right\} d\Phi \left( c_{i}| \Sigma \right) \,,
\label{CH12X_eq_16} \\
\text{with }L_{i}\left( \Omega _{-\Sigma }| c_{i}\right)
&=&\prod_{m=1}^{M}\prod_{t=1}^{T}\prod_{j=1}^{J}\left[ P_{imt}R_{imtj}\right]
^{D_{imt}E_{ijt}}\,,  \notag
\end{eqnarray}%
where $\Phi \left( .\right) $ is the multi-dimensional cumulative normal
distribution and $L_{i}\left( \Omega _{-\Sigma }| c_{i}\right) $ is the
likelihood function, for household $i$ and for a given value of $\Omega
_{-\Sigma }$ and $c_{i}$. $D_{imt}$ (resp. $E_{ijt}$) is a dummy variable
taking a value of 1 if household $i$ chooses to move for reason $m$ (resp.
chooses to stay for duration category $j$) on the $t$-th choice occasion and
0 otherwise. Equation (\ref{CH12X_eq_16}) needs the evaluation of a
multi-dimensional integral of size equal to the number of rows in $c_{i}$.
\citeA{Eluru2009} apply \index{Quasi-Monte Carlo}Quasi-Monte Carlo simulation techniques
based on the \index{Halton sequence}Halton sequence to approximate this integral in the likelihood
function and maximize the logarithm of the resulting simulated likelihood
function across individuals with respect to $\Omega $ (see \citeA{Bhat2001,
Bhat2003}). \citeA{Eluru2009} use a longitudinal data set of households from
a stratified sample of municipalities in the \index{housing!model example!Zurich}Zurich region of Switzerland
over the period 1985-2004. The data set includes 1012 households and 2590
move records. They found that several demographic, socioeconomic, and
commute related variables (\textit{e.g.}, age, gender, family reasons,
education/employment reasons, accommodation related reasons, surrounding
environment related reasons, vicinity to family and friends, etc.) have a
significant influence on the reason for move and the \index{duration!of stay}duration of stay. In
the \index{duration!of stay}duration of stay model, \citeA{Eluru2009} found that household size
creates \index{heterogeneity}heterogeneity across the sample of households. They show that people
who own \index{dwelling}dwellings have a lower probability of moving for surrounding
vicinity related reasons than those renting their units. Likewise, people
who live in smaller homes have higher probabilities of short duration stays
probably because they are looking for larger homes. Having a mix of job
opportunities located close to residential \index{neighborhood}neighborhoods increases the
\index{duration!of stay}duration of stay in the \index{dwelling}dwelling. Reducing commute distances promotes longer
durations of stay, etc. \citeA{Eluru2009} found that common unobserved
factors jointly affect the reason to move and the \index{duration!of stay}duration of stay and calls
for a joint modeling framework that allows error correlation structures.

\index{endogeneity}Endogeneity (or simultaneity) is a fundamental aspect of modelling housing
that should be taken into account both for hedonic \index{housing!price}housing price functions
and for choice models of \index{residential!location}residential location. This is the object of the
next section.

\section{Multi-dimensional Dynamic Models of Housing Models}

\label{CH12X_section5}

In hedonic \index{housing!price}housing price functions, some explanatory variables, in addition
to the dependent variable and its \index{spatial!lag}spatial lag, may be endogenous
following the simultaneous choice of the \index{housing!price}house price and of the quantities
of \index{attributes}attributes. This is particularly true for floor space (see \citeA
{Fingleton2008b} who extended \citeA{Kelejian1998} feasible \index{generalized!spatial two-stage least squares}generalized spatial two-stage least squares estimator to account for
endogenous variables due to system feedback, given an autoregressive or a
moving average error process). As for \index{hedonic!price}hedonic price functions, \index{endogeneity}endogeneity
is expected to occur mainly as a result of the omission of \index{attributes}attributes in
\index{discrete choice}discrete choice models of \index{residential!mobility}residential mobility. In the literature, several
methods have been proposed to consider \index{endogeneity}endogeneity. \citeA{Berry1995}
proposed a \index{fixed effects}fixed effects procedure, by product and market, to solve
market-level \index{endogeneity}endogeneity in the automobile sector. \citeA{Guevara2006}
applied to residential \index{location choice}location choice models the \index{control function method}control function method,
which is based on the inclusion of an additional variable that controls for
the \index{endogeneity}endogeneity problem (see \citeA{Heckman1978} and \citeA{Blundell2004}).
They applied residential \index{location choice}location choice models based on 630 households of
renters who had moved to their present location between 1999 and 2001 in
\index{housing!model example!Santiago (Chile)}Santiago (Chile). The results show that price \index{endogeneity}endogeneity is significant in
choice models of \index{residential!location}residential location and that the \index{control function method}control function method
is suitable to account for it.

\index{endogeneity}Endogeneity is not limited to the correlation between the dependent
variables and \index{attributes}attributes (in the equation or omitted) or to the simultaneity
of demand and supply, the marginal \index{willingness!to pay}willingness to pay and the marginal
\index{willingness!to accept}willingness to accept. \index{location choice}Location choices and \index{housing!investments}housing investments are
fundamentally dynamic decisions over multiple time periods. In the 2D panel
data literature, some \index{dynamic model}dynamic models have been applied on \index{real estate}real estate
topics. For instance, \citeA{Engle1985} used a version of a dynamic \index{multiple-indicator multiple-cause}multiple-indicator multiple-cause (DYMIMIC) model
for an \index{hedonic!price}hedonic price
model of the \index{resale housing market}resale housing market for a suburb of \index{housing!model example!US San Diego}San Diego, California
during the period 1973-1980. The specification of the model features hedonic
equations for each house sale and a dynamic equation for the capitalization
rate which is taken to be an unobservable time series to be estimated
jointly with the unknown parameters. \citeA{Engle1985} used \index{maximum likelihood}maximum likelihood with an \index{EM algorithm}EM algorithm
based upon \index{Kalman filtering}Kalman filtering. Some authors
have used, in a 2D framework, the \index{dynamic factor models}dynamic factor models (DFM)
and/or \index{large-scale Bayesian vector autoregressive}{large-scale Bayesian vector autoregressive (LBVAR) models to
forecast \index{housing!price}housing price. These models are interesting to study the
\index{ripple effect}\textquotedblleft ripple effect", \textit{i.e.}, the propagation of
shocks to \index{housing!price}house prices across regions. For instance, \citeA{Das2010} forecast
regional \index{housing!price}house price inflation for five metropolitan areas of \index{housing!model example!South Africa}South Africa
using principal components obtained from quarterly macroeconomic time series
in the period 1980 to 2006. In the majority of the cases, the dynamic factor
model statistically outperforms the vector autoregressive models, using both
the classical and the \index{Bayesian!treatments}Bayesian treatments. They also considered spatial and
non-spatial specifications. \citeA{Das2010} indicate that macroeconomic
fundamentals in forecasting \index{housing!price}house price inflation are important. \citeA
{Li2011} investigate moving trends of \index{housing!price}house prices in 42 metropolitan areas
in the \index{housing!model example!US}United States from the perspective of large-scale models, which are
also DFM and LBVAR models. These models accommodate a large panel data
comprising 183 monthly series for the U.S. economy, and an in-sample period
of 1980 to 2007 are used to forecast one- to twelve-months-ahead \index{housing!price}house price
growth rate over the out-of-sample horizon of 2008 to 2010. \citeA{Li2011}
show that \index{dynamic factor models}DFM consistently outperforms its \index{large-scale Bayesian vector autoregressive}LBVAR alternative for forecasting
the \index{housing!price}house price growth rate for the overall U.S. housing market. The
forecasting power of \index{dynamic factor models}DFM does not decrease as the number of forecasted
period ahead increases, while \index{large-scale Bayesian vector autoregressive}LBVAR has its best performance for
two-month-ahead forecast and then its forecasting accuracy decays. \citeA
{Beenstock2015} using data from 9 regions of \index{housing!model example!Israel}Israel over 1987-2010 apply
spatial \index{panel cointegration}panel cointegration methods for a \index{dynamic model}dynamic model of regional housing
markets in which people prefer to live where housing is cheaper and building
contractors prefer to build in regions where construction is more
profitable. Based on dynamic \index{hedonic!price}hedonic price functions, the analysis of
nonstationary spatial panel data shows that although housing starts vary
directly with profitability as measured by \index{housing!price}house prices relative to building
costs, they vary inversely with profitability in neighboring regions. \citeA
{Beenstock2015} show that there is a non negligible spatial substitution in
housing construction and this substitution effect suggests that contractors
have local building preferences since they regard neighboring regions as
close substitutes but not more distant regions. \citeA{Abate2016} investigate
the interactions between \index{housing!price}house price fluctuations and output growth rate
across 373 metropolitan statistical areas in the US over the period
2001-2013. In a panel data context, they use time varying spatial
econometric \index{hedonic!price}hedonic price functions. They show that the spatial correlation
coefficient across metropolitan areas has been increasing over time,
indicating an increasing synchronization of \index{housing!price}house prices across metropolitan
statistical areas during the sample period.

\index{spatio-temporal}Spatio-temporal models of \index{hedonic!price}hedonic price functions have been
recently proposed to jointly take into account time effects and spatial
effects either through multifactor error structure or through specific
weight matrices. For instance, \citeA{Holly2010} considered the determination
of real \index{housing!price}house prices in a panel made up of 49 \index{housing!model example!US}US States over 29 years. An
error correction model with a cointegrating relationship between
real \index{housing!price}house prices and real incomes is found once they take proper account of
both \index{heterogeneity}heterogeneity and cross sectional dependence. See also \citeA{Latif2015}
for a study on the impact of new immigration on housing rent, using \index{housing!model example!Canada}Canadian
province-level panel data from 1983 to 2010. \citeA{Latif2015} uses \index{panel cointegration}panel cointegration regressions and panel vector error correction models
and shows that immigration flow has a significant positive impact on housing
rent both in the short run and in the long run. There are also extensions of
the spatial \index{hedonic!price}hedonic price functions which use a weight matrix that expresses
\index{spatio-temporal}spatio-temporal relations instead of purely spatial. A general $(N\times N)$
\index{spatio-temporal}spatio-temporal weight matrix $W$ is obtained by splitting its construction
into two separate matrices of the same dimension. The first matrix, $S$,
captures the spatial relations among the $N$ observations and a second
matrix, $T$, expresses the temporal direction of observations. \citeA
{Smith2011} have proposed a \index{spatio-temporal}spatio-temporal weight matrix defined as the
Hadamard product between two spatial and temporal distance weight matrices $%
W=S\odot T=\left[ s_{jl}\right] \odot \left[ t_{jl}\right] $. It identifies
the \index{spatio-temporal}spatio-temporal \index{neighbor}neighbors that affect \index{hedonic!price}hedonic price determination. The
elements $s_{jl}$ indicate the way observation $j$ is spatially connected to
observation $l$. The elements on the diagonal $s_{jj}$ are set to zero while
the off-diagonal elements are defined by an inverse distance function: $%
s_{jl}=d_{jl}^{-\gamma }$ if $d_{jl}<\overline{d}$ and 0 elsewhere where $%
d_{jl}$ is the geographic distance between locations $j$ and $l$, $d_{jl}<%
\overline{d}$ is a critical cut-off and $\gamma \geq 0$. The elements $t_{jl}
$ represent the time elapsed between realization of observations $j$ and $l$%
. One assumes that observations have been ordered chronologically so the
first row of $T$ corresponds to the earliest observation while the last row
corresponds to the latest observation. The elements on the diagonal $t_{jj}$
are set to zero while the off-diagonal elements are defined by an inverse
function of the time elapsed between two observations: $t_{jl}=|
t_{j}-t_{l}| ^{-\alpha }$ if $| t_{j}-t_{l}| <\overline{t}$ and 1
elsewhere. $t_{j}$ (resp. $t_{l}$) is the time when \index{dwelling}dwelling $j$ (resp. $l$)
is sold. $\overline{t}$ is a critical cut-off value and $\alpha $ is a
penalty parameter to be fixed. Several authors have used \index{spatio-temporal}spatio-temporal
models of \index{hedonic!price}hedonic price functions with standard spatial specifications (%
\index{spatial!autoregressive}spatial autoregressive (SAR), \index{spatial!error}spatial error (SEM), \index{spatial!Durbin model}spatial Durbin model, etc.)
but with different \index{spatio-temporal}spatio-temporal matrices $W$. They got
better results in terms of estimation and/or forecasting as compared to
those obtained with usual purely \index{spatial!weight}spatial weight matrices. See for instance,
\citeA{Pace2000} for an application on the residential market of \index{housing!model example!US B\^{a}ton Rouge}B\^{a}ton Rouge, Louisiana, during 1984-1992, \citeA{Liu2013} for an application of
housing in \index{housing!model example!The Netherlands}Randstad, The Netherlands, during the years 1997-2007, \citeA
{Nappi-Choulet2011} for the residential market of \index{housing!model example!Paris}Paris for the years
1995-2005, or \citeA{Thanos2016} for the \index{housing!model example!Aberdeen (Scotland)}Aberdeen, Scotland, housing market
during 2004-2007, to mention a few. Unfortunately, and to our knowledge,
nobody has used these \index{spatio-temporal}spatio-temporal multifactor error structures or the
\index{spatio-temporal}spatio-temporal weight matrices in a \index{three-dimensional}three-dimensional framework. But, it
could be a promising development for future research.

The developments in the dynamics of modelling housing is not only focused on
\index{hedonic!price}hedonic price functions. Some authors have been interested in dynamic
versions of discrete models of \index{location choice}location choice. \index{forward-looking}Forward-looking
behavior in the housing market justify dynamic considerations in a model of
\index{location choice}location choice. Several authors have underlined the need to use dynamic
specifications of modelling housing. For instance, \citeA{Case2012}, using
questionnaire surveys for home buyers in four \index{housing!model example!US}U.S. cities over 2003-2012,
have shown that the root causes of the speculative bubble can be seen in
their long-term home price expectations, which reached abnormal levels
relative to the \index{mortages}mortgage rate at the peak of the boom and declined sharply
since. The downward turning point around 2005 of the long boom that preceded
the crisis was associated with changing public understanding of speculative
bubbles. But estimating dynamic discrete models of \index{location choice}location choice is a
rather challenging and stimulating objective and is technically difficult.
\citeA{Bayer2016} noted that first, estimation of residential \index{sorting}sorting and
\index{hedonic!equilibrium model}hedonic equilibrium models needs to match a large sample of households,
their characteristics to the location and the features of their housing
choices. Second, the high dimensionality of the state space (consisting of
current lifetime utilities and \index{neighborhood}neighborhood characteristics) -- required to
define the evolution of a urban system -- leads to the curse of
dimensionality which puts a brake on the estimation of an acceptable sized
\index{dynamic model}dynamic model of \index{residential!location}residential location decisions.

\citeA{Diao2015} propose a \index{real option theory}real-option based \index{dynamic model}dynamic model to simulate real
estate developer behavior. In a \index{three-dimensional}three-dimensional framework (property, type
of property and time for private residential housing in \index{housing!model example!Singapore}Singapore during
1995-2012), they extend the standard \index{discrete choice}discrete choice model approach by
adding an explicit probabilistic representation of development templates
available to developers to take into account both developers' option to hold
the land undeveloped and the market volatility of different development
types. In their proposed simulation framework, \citeA{Diao2015} suppose that
a developer making investment decisions for a parcel faces a set of
alternative development templates in a market with uncertainty. At each time
period, the developer estimates future revenue and construction cost of
feasible development templates under planning constraints and related real
option values. He chooses the template based on the principle of profit
maximization, but only does so if the return of the development template is
greater than a threshold level (value of the call option), which is a
function of the market volatility of the built property as suggested by the
\index{real option theory}real option theory, otherwise, keeps the \textit{status quo}. The model
components in the proposed simulation framework are calibrated with private
housing data in \index{housing!model example!Singapore}Singapore. The results show significant volatility in
\index{housing!price}housing price and construction costs, relevant differences in volatility
across housing types, and good fit, in the hedonic model, of market prices
and construction costs. This kind of research contributes to the
microsimulation literature by proposing an interesting approach which takes
into account the dynamic and volatile nature of the \index{real estate}real estate market but,
unfortunately, this remains a simulation study.

\citeA{Bayer2016} have proposed a new approach for estimating a
\index{three-dimensional}three-dimensional \index{dynamic model}dynamic model of demand for houses and \index{neighborhood}neighborhoods that
is computationally tractable. Using a \index{semi-parametric estimation}semi-parametric estimation
approach, they control for unobserved household and \index{neighborhood}neighborhood
\index{heterogeneity}heterogeneity. Their model adapts dynamic demand models for durable goods in
a housing market context. They treat houses as assets and allow households
wealth to evolve endogenously. Households anticipate selling their homes at
some point in the future and then consider the expected evolution of
\index{neighborhood}neighborhood \index{amenities}amenities and \index{housing!price}housing price when deciding where and when to
purchase, or sell, their house. They relax the \index{index sufficiency assumption}index sufficiency assumption
which is standard in the dynamic demand literature. This assumption helps to
deal with the computational challenges posed by the large state space
typically arising in models of dynamic demand. Instead of treating the logit
\index{inclusive value}inclusive value as a sufficient statistic for predicting future continuation
values, \citeA{Bayer2016} define the continuation value from predicted future
lifetime utilities, which depend on the state space in a flexible manner.
Last, they use stable and uniform realtor fees to estimate the marginal
utility of consumption without the need for a price instrument. They use the
fact that households face a monetary trade-off both in the standard sense of
deciding which product (\index{neighborhood}neighborhood) to purchase but also in terms of
deciding when to move. They take advantage of the fact that realtor fees
during the sample period were quite uniform ($6\%$ of the house value) in
order to identify the marginal utility of consumption when estimating each
resident's move-stay decision. The decision variable, $d_{it}$, denotes both
of the choices made by household $i$ in period $t$, whether to move and
where to move, conditional on deciding to move. If a household decides to
move, the decision is denoted $d_{it}=j$, $j=0,1,...,J$ where $j$ indexes
\index{neighborhood}neighborhoods, $J$ denotes the total number of \index{neighborhood}neighborhoods in the region
and $0$ denotes the outside option. The data concern housing transactions in
the \index{housing!model example!US San Francisco}San Francisco Bay Area from 1994-2004 for more than $220,000$ households
and $2398$ \index{neighborhood}neighborhoods. We give only some results as the paper is highly
technical. But, the model and estimation procedure presented in this paper
are very general and can be applied to a broad range of dynamic studies in
housing markets. The model uses a two-stage estimator. In the first stage,
\citeA{Bayer2016} use the household location and the mobility decisions to
estimate the value of lifetime expected utility for each \index{neighborhood}neighborhood, time
period, and household type as well as an unobservable characteristic that
captures a household's preference for sub-regions within the \index{housing!model example!US San Francisco}San Francisco
Bay Area. In the second stage, they recover fully-flexible estimates of
per-period utility and regress them on a set of observable \index{attributes}attributes. They
use a \index{semi-parametric estimation}semi-parametric estimation approach to control for the \index{endogeneity}endogeneity of
price in this second stage, utilizing outside information relating to the
financial cost of moving to pin down the coefficient on \index{housing!price}house prices. The
results indicate that the downward biases associated with static demand
estimation are significant for three important non-marketed \index{amenities}amenities: air
\index{quality}quality, crime, and \index{neighborhood}neighborhood race. For instance, for a $10\%$ change in
each amenity, the static model overestimates the \index{willingness!to pay}willingness to pay for
living in close proximity to \index{neighbor}neighbors of the same race for low-income
households. The static estimation is $\$1,627.03$ whereas the corresponding
dynamic estimation is $\$612.09$. For high-income households, the bias runs
in the opposite direction and the static model underestimates the
\index{willingness!to pay}willingness to pay by a factor of more than two. The static model always
underestimates the \index{willingness!to pay}willingness to pay for living in close proximity to
crimes. For low-income households and for a $10\%$ increase in violent
crime, the static estimation is -$\$291.14$ while the corresponding dynamic
estimation is -$\$350.18$. This is also true for air pollution.


\section{From \index{variational inference}Variational Inference ... to Multi-dimensional Housing Models}

\label{CH12X_section6}

The presence of numerous \index{latent variables}latent variables, omitted variables, the definition of dynamic and spatial
structures within multi-dimensional frameworks (3D, 4D or more) and the
econometric complexity that results will not make things easier and must
move us towards the use of flexible models and methods. Among many
others promising future pathways is probably the use of \index{variational inference!variational Bayesian approximations}variational Bayesian approximations
 (see for instance \citeA{Ormerod2010}, \citeA{Lee2016}, \citeA{blei2017}, \citeA{baltagi2019carbon}, \citeA{baltagi2020growth} or \citeA{nolanwand2020} to mention a few).
These methods facilitate approximate inference for the parameters in complex
statistical models and provide fast, deterministic alternatives to Monte
Carlo methods to potentially overcome many problems in applied modelling of
housing. While the two-level approaches have been developed over the past decade, the three-level approach was recently proposed by \citeA{nolanwand2020}.
Their \index{variational inference!streamlined variational inference}streamlined variational inference algorithm for \index{variational inference!three-level random effects model}three-level random effects models relies on theorems provided by \citeA{nolan2020} concerning linear system solutions and sub-blocks of matrix inverses for  \index{variational inference!three-level sparse matrix least squares problem}three-level sparse matrix problems which are the basis for their  \index{variational inference!QR-decomposition-based streamlined algorithm}QR-decomposition-based streamlined algorithm. This algorithm allows one to obtain \index{variational inference!mean field variational Bayes approximation (MFVB)}mean field variational Bayes approximate posterior density functions for the parameters in the \index{variational inference!three-level linear mixed model}three-level linear mixed model. We will briefly present their method and apply it to data on housing transaction prices in \index{housing!model example!Paris}Paris in 2003 for 20 {\it arrondissements}, each divided in four {\it quartiers}. Even though four-level and even higher level situations may be of interest, unfortunately the required theory is not yet in place.\\
Inference based on MCMC can be very slow for
such  two-level or three-level linear mixed model and MCMC methods may suffer from poor mixing. Variational Bayesian inference can help in tackling the scalability
challenge of big data sets and/or models with large sparse covariance matrices as they use a deterministic optimization approach to approximate
the posterior distribution. The parameters of the approximate distribution are chosen to minimize some measure of distance (as the \index{variational inference!Kullback-Leibler divergence}Kullback-Leibler
divergence) between the approximation and the posterior. \index{variational inference!mean field variational Bayes approximation (MFVB)}Mean field
variational Bayes approximation is analogous to Gibbs sampling for conjugate models (see \citeA{bishop2006}, \citeA{Ormerod2010}, \citeA{pham2013} and \citeA{Lee2016} to mention a few). \\
Consider a generic Bayesian model with observed vector $y$ and parameter vector $\theta$ that is continuous over the parameter space $\Theta$. The Bayes theorem allows one to define the posterior distribution
as: 
\begin{equation*}
p\left( \theta | y\right) =\frac{p\left( \theta ,y\right) }{p\left(
y\right) }=\frac{p\left( y | \theta \right) p\left( \theta \right) }{%
p\left( y\right) }\text{ with }p\left( y\right) =\int_{\Theta }p\left(
\theta ,y\right) d\theta  \label{CH12X_eq_17}
\end{equation*}%
Let $q$ be an arbitrary density function over $\Theta $. Then, the
logarithm of the marginal likelihood satisfies (see \citeA{bishop2006}, \citeA{Ormerod2010}):
\begin{eqnarray*}
\log p\left( y\right) &=&\log p\left( y\right) \int_{\Theta }q\left( \theta
\right) d\theta =\int_{\Theta }q\left( \theta \right) \log p\left( y\right)
d\theta  \notag \\
&=&\int_{\Theta }q\left( \theta \right) \log \left\{ \frac{p\left( \theta
,y\right) /q\left( \theta \right) }{p\left( \theta | y\right) /q\left(
\theta \right) }\right\} d\theta  \notag \\
&=&\int_{\Theta }q\left( \theta \right) \log \left\{ \frac{p\left( \theta
,y\right) }{q\left( \theta \right) }\right\} d\theta +\int_{\Theta }q\left(
\theta \right) \log \left\{ \frac{q\left( \theta \right) }{p\left( \theta
| y\right) }\right\} d\theta  \notag \\
&=&\log \underline{p}\left( y,q\right) +KL(q,p)  \label{CH12X_eq_18}
\end{eqnarray*}
where $KL(q,p)$ is the Kullback-Leibler divergence between $q\left( \theta
\right) $ and $p\left( \theta | y\right) $. Furthermore, $\log \underline{%
p}\left( y,q\right) $ is a lower bound on the marginal log-likelihood. The
Kullback-Leibler divergence becomes 
\begin{eqnarray*}
KL(q,p) &=&E_{q}\left[ \log q\left( \theta \right) %
\right] -E_{q}\left[ \log p\left( \theta | y\right) %
\right]  \notag \\
&=&E_{q}\left[ \log q\left( \theta \right) \right]
-E_{q}\left[ \log p\left( \theta ,y\right) \right]
+\log p\left( y\right)   \label{CH12X_eq_19}
\end{eqnarray*}%
where the last term, $\log p\left( y\right) $, is a constant. The
minimization of the Kullback-Leibler divergence is thus equivalent to
maximizing the scalar quantity, 
\begin{equation*}
\log \underline{p}\left( y,q\right) =E_{q }\left[ \log
\left( \frac{p\left( \theta ,y\right) }{q\left( \theta \right) }\right) %
\right]   \label{CH12X_eq_20}
\end{equation*}%
which is usually referred as the \index{variational inference!evidence lower bound (ELBO)}evidence lower bound (ELBO).\\
 Let $\left\{ \theta _{1},...,\theta _{M}\right\} $ be a partition of
the parameter vector $\theta $. The \index{variational inference!mean field variational Bayes approximation (MFVB)}MFVB approximates the posterior
distribution $p\left( \theta | y\right) $ by the product of the \index{variational inference!$q$-densities}$q$%
-densities:\footnote{
This is known as the \textit{mean field restriction}. The term \textit{mean
field} originated from physics.} 
\begin{equation*}
q\left( \theta \right) =\prod_{j=1}^{M}q_{j}\left( \theta _{j}\right)
 \label{CH12X_eq_21a}
\end{equation*}%
The \index{variational inference!optimal $q$-densities}optimal $q$-densities which minimize the
Kullback-Leibler divergence are given by 
\begin{equation}
q_{j}^{\ast }\left( \theta _{j}\right) \propto \exp \left[ E_{q\left(
-\theta _{j}\right) }\left\{ \log p\left( \theta _{j}|\text{ rest}\right)
\right\} \right] \text{ , }j=1,...,M   \label{CH12X_eq_21}
\end{equation}%
where $E_{q\left( -\theta _{j}\right) }$ denotes expectation with respect to 
$\prod_{k\neq j}q_{k}\left( \theta _{k}\right) $.\newline
$\text{rest}\equiv \left\{ y,\theta _{1},...,\theta _{j-1},\theta
_{j+1},...,\theta _{M}\right\} $ is the set containing the rest of the
random vectors in the model, except $\theta _{j}$ and the distributions $%
\left( \theta _{j}|\text{ rest}\right) $ are the full conditionals in the
MCMC literature.\\
The iterative scheme for obtaining the \index{variational inference!optimal $q$-densities}optimal $q$-densities under
product restriction (\ref{CH12X_eq_21}) is\\ 
\noindent \rule{12cm}{1pt}
\begin{enumerate}
\item Initialize $q^{*}_1 \left( \theta _{1}\right) \text{ , } q^{*}_2
\left( \theta _{2}\right) \text{ , } \cdots \text{ , } q^{*}_M \left( \theta
_{M}\right)$
\item Cycle through updates: 
\begin{eqnarray*}
q^{*}_1 \left( \theta _{1}\right) & \leftarrow & \frac{ \exp \left[
E_{q\left( -\theta _{1}\right) }\left\{ \log p\left( y , \theta \right)
\right\} \right] }{ \int \exp \left[ E_{q\left( -\theta _{1}\right) }\left\{
\log p\left( y , \theta \right) \right\} \right] d \theta_1 } \\
\cdots & \cdots & \cdots \cdots \cdots \cdots \cdots \cdots \\
q^{*}_M \left( \theta _{M}\right) & \leftarrow & \frac{ \exp \left[
E_{q\left( -\theta _{M}\right) }\left\{ \log p\left( y , \theta \right)
\right\} \right] }{ \int \exp \left[ E_{q\left( -\theta _{M}\right) }\left\{
\log p\left( y , \theta \right) \right\} \right] d \theta_M }
\end{eqnarray*}
until the increase in $\log \underline{p}\left( y,q\right)$ is negligible.
\end{enumerate}
\noindent \rule{12cm}{1pt}

Compared to the minimization of the KL divergence, the maximization of the
\index{variational inference!evidence lower bound (ELBO)}ELBO is often a more convenient objective of the optimization over the free
distributional parameters. \citeA{nolanwand2020} apply this principle and derive the \index{variational inference!mean field variational Bayes approximation (MFVB)}MFVB
approximation of the \index{variational inference!three-level linear mixed model}three-level linear mixed model.\\
Suppose that the response data vector $y$ is modeled according to a Bayesian version of the \index{variational inference!Gaussian linear mixed model}Gaussian linear mixed model\footnote{This draws heavily from the results in \citeA{nolanwand2020}.} 
\begin{eqnarray*}
y | \beta, u, R \sim N\left( X \beta  + Z u , R    \right) \text{, } u | G \sim N(0,G)  \text{, } \beta \sim N( \mu_{\beta} , \Sigma_{\beta}  )  \label{CH12X_eq_22}
\end{eqnarray*}
for hyperparameters $\mu_{\beta}$ and $\Sigma_{\beta}$ and such that $\beta$ and $u | G$ are independent. $X$ and $Z$ are the fixed effects and random effects design matrices
associated with the fixed effects and random effects vectors $\beta $ and $u$.\footnote{This terminology is different from what the panel data literature dubs as
\textquotedblleft fixed\textquotedblright\ and \textquotedblleft random\textquotedblright\ effects.} Prior specification for the covariance matrices $G$ and $R$ involves auxiliary covariance matrices $A_G$ and $A_R$ with conjugate \index{variational inference!Inverse G-Wishart distributions}Inverse G-Wishart distributions (generalization of inverse Wishart distributions)\footnote{See Appendix.} (\citeA{wand2017fast}).\\ 
Full Bayesian inference for the $\beta$, $G$ and $R$ and the random effects $u$ involves the posterior density function $p(\beta, u, A_G, A_R, G,R | y)$, but typically is analytically intractable and Markov chain Monte Carlo approaches are required for practical exact inference.  Variational approximate inference involves mean field restrictions such as
\begin{eqnarray}
p\left( \beta, u, A_G, A_R, G,R | y \right) & \approx & q \left( \beta, u, A_G, A_R \right) q\left( G,R \right) \notag \\
     & \approx &  q \left( \beta, u \right)  q( A_G) q(A_R) q( G) q(R)  \label{CH12X_eq_23}
\end{eqnarray}
Then, the forms and optimal parameters for the $q$-densities are obtained by minimizing the \index{variational inference!Kullback-Leibler divergence}Kullback-Leibler divergence of the right-hand side of (\ref{CH12X_eq_23}) from its left-hand side. The \index{variational inference!optimal $q$-densities}optimal $q$-density parameters are interdependent and a coordinate ascent algorithm (see \citeA{Ormerod2010}) is used to obtain their solution. The optimal $q$-density for $(\beta,u)$, denoted by $q^{*}(\beta,u)$, is a multivariate Normal density function with mean vector $\mu_{q(\beta,u)}$ and covariance matrix $\Sigma_{q(\beta,u)} $. The coordinate ascent algorithm is such that they are updated according to (see \citeA{Lee2016}, \citeA{nolanwand2020})\footnote{Very often, $R=\sigma^{2}_{\varepsilon} I$ where $I$ is an identity matrix and where $(\sigma^{2}_{\varepsilon} | a_{\sigma^{2}_{\varepsilon}})$ follows an inverse $\chi^{2}$ and $a_{\sigma^{2}_{\varepsilon}} (\equiv A_R)$ also follows an inverse $\chi^{2}$ distribution.}
\begin{eqnarray*}
\Sigma_{q(\beta,u)} & \leftarrow & \left\{  C^{\prime} E_{q}(R^{-1})  C + 
\left[
\begin{array}{cc}
\Sigma^{-1}_{\beta} & 0 \\
0    & E_{q}(G^{-1})
\end{array}
\right]
  \right\}    \notag \\
 \mu_{q(\beta,u)} & \leftarrow & \Sigma_{q(\beta,u)} C^{\prime} E_{q}(R^{-1}) \left( y + 
 \left[
\begin{array}{c}
\Sigma^{-1}_{\beta} \mu_{\beta} \\
0    
\end{array}
\right]
 \right)
 \label{CH12X_eq_24}
\end{eqnarray*}
where $E_{q}(G^{-1})$ and $E_{q}(R^{-1})$ are the \index{variational inference!$q$-density expectations}$q$-density expectations of $G^{-1}$ and $R^{-1}$ and $C=[X, Z]$.\\
Let us now turn our attention to the \index{variational inference!mean field variational Bayes approximation (MFVB)}mean field variational Bayes (MFVB) approximation of the \index{variational inference!three-level linear mixed model}three-level linear mixed model proposed by \citeA{nolanwand2020}.
\begin{eqnarray*}
y_{ij} | \beta, u^{L1}_{i}, u^{L2}_{ij}, \sigma^{2} \sim N\left( X_{ij} \beta  + Z^{L1}_{ij} u^{L1}_{i} + Z^{L2}_{ij} u^{L2}_{ij}, \sigma^{2}_{\varepsilon} I    \right) 
 \label{CH12X_eq_25}
\end{eqnarray*}
where $1\leq i \leq m$ is the index of the outer group (the 20 {\it arrondissements} of Paris) and $1\leq j \leq n_i$ is the index of the inner group (the 4 {\it quartiers} of each {\it arrondissement} of Paris in our particular example). Inside each inner group, we have $o_{ij}$ observations for each variable ($o_{ij}$ housing transaction prices in our example). Of course, we can have unbalanced inner groups with, in each, unbalanced observations for each variable.  $y_{ij}$ is a $\left( o_{ij} \times 1\right)$ vector. $X_{ij}$ is an $(o_{ij}  \times p)$ matrix of covariates, $Z^{L1}_{ij}$ (resp.  $Z^{L2}_{ij}$) is an $(o_{ij} \times q_1)$ (resp. $(o_{ij} \times q_2)$) block-diagonal matrix of the $X_{ij}$ submatrices.\footnote{If $q_1\neq q_2$, it means that some random effects of level 1 or level 2 concern only some covariables. If $p \neq q_1$ (with $q_1 < p$) or $p \neq q_2$ (with $q_2 < p$), it also means that some covariables have only fixed effects. In our example, we will assume that $p=q_1=q_2$.}
Thus, we have a \index{variational inference!three-level linear mixed model}three level random-intercept-random-slopes model with the first level (outer group), the second level (inner group) and the third level being the set of observations in each inner group.\\
The full Bayesian model (with priors on parameters and
hyperparameters) is given by (see \citeA{nolanwand2020}):
\begin{eqnarray*}
y_{ij} & | & \beta, u^{L1}_{i}, u^{L2}_{ij}, \sigma^{2} \sim N\left( X_{ij} \beta  + Z^{L1}_{ij} u^{L1}_{i} + Z^{L2}_{ij} u^{L2}_{ij}, \sigma^{2}_{\varepsilon} I    \right)   \label{CH12X_eq_26} \\
&&\text{with}  \notag \\
&&\left\{ 
\begin{array}{lcl}
\left(
\begin{array}{c}
u^{L1}_{i} \\
u^{L2}_{ij}
\end{array}
\right)
| \Sigma ^{L1}, \Sigma ^{L2}   \sim  N\left( 
\left(
\begin{array}{c}
0 \\
0
\end{array}
\right) ,
\left( 
\begin{array}{cc}
\Sigma ^{L1} & 0 \\
0   & \Sigma ^{L2}
\end{array}%
\right) \right) \\ 
\beta \sim  N\left( \mu_{\beta},\Sigma _{\beta } \right) \text{, }
\sigma^{2}_{\varepsilon} | a_{\sigma^{2}_{\varepsilon}} \sim \text{Inverse-} \chi^{2} (\nu_{\sigma^{2}_{\varepsilon}}, 1/a_{\sigma^{2}_{\varepsilon}} ) \\ 
a_{\sigma^{2}_{\varepsilon}} \sim \text{Inverse-} \chi^{2} (1, 1/( \nu_{\sigma^{2}_{\varepsilon}} s^{2}_{\sigma^{2}_{\varepsilon}}   )) \\ 
\Sigma ^{L1}| A_{\Sigma ^{L1}} \sim \text{Inverse-G-Wishart}(G_{full}, \nu_{\Sigma ^{L1}} + 2q_1 -2, ( A_{\Sigma ^{L1}})^{-1}   )\\ 
A_{\Sigma ^{L1}} \sim  \sim \text{Inverse-G-Wishart}(G_{diag}, 1, \left\{  \nu_{\Sigma ^{L1}} diag( s^{2}_{\Sigma ^{L1},1}, \cdots ,  s^{2}_{\Sigma ^{L1},q_1}   )   \right\}^{-1}     )\\ 
\Sigma ^{L2}| A_{\Sigma ^{L2}} \sim \text{Inverse-G-Wishart}(G_{full}, \nu_{\Sigma ^{L2}} + 2q_2 -2, ( A_{\Sigma ^{L2}})^{-1}   )\\ 
A_{\Sigma ^{L2}} \sim  \sim \text{Inverse-G-Wishart}(G_{diag}, 1, \left\{  \nu_{\Sigma ^{L2}} diag( s^{2}_{\Sigma ^{L2},1}, \cdots ,  s^{2}_{\Sigma ^{L2},q_2}   )   \right\}^{-1}     )\
\end{array}
\right.  \notag
\end{eqnarray*}
and $ \nu_{\Sigma ^{L1}}>0$, $\nu_{\Sigma ^{L2}}>0$, $s^{2}_{\Sigma ^{L1},1}, \cdots ,  s^{2}_{\Sigma ^{L1},q_1} > 0$ and $ s^{2}_{\Sigma ^{L2},1}, \cdots ,  s^{2}_{\Sigma ^{L2},q_2} > 0$. The minimal mean field restriction needed for a tractable variational inference algorithm is 
\begin{eqnarray*}
p( \beta , u,   a_{\sigma^{2}_{\varepsilon}},  A_{\Sigma ^{L1}}, A_{\Sigma ^{L2}} , \sigma^{2}_{\varepsilon} , \Sigma ^{L1}, \Sigma ^{L2} | y  ) \approx
q(\beta , u,   a_{\sigma^{2}_{\varepsilon}},  A_{\Sigma ^{L1}}, A_{\Sigma ^{L2}}) q(\sigma^{2}_{\varepsilon} , \Sigma ^{L1}, \Sigma ^{L2}  ) \notag \\
 \label{CH12X_eq_27}
\end{eqnarray*}
The \index{variational inference!optimal $q$-densities}optimal $q$-density functions for the parameters of interest being as follows:
\begin{equation*}
\left\{
\begin{array}{l}
q^{*}(\beta, u) \text{ has a } N(\mu_{q(\beta , u)}, \Sigma_{q(\beta , u)} ) \text{ density function} \\
q^{*}( \sigma^{2}_{\varepsilon})  \text{ has an Inverse-} \chi^{2} (\xi_{q(\sigma^{2}_{\varepsilon})}, \lambda {q(\sigma^{2}_{\varepsilon})} ) \text{ density function}  \\
q^{*}(  \Sigma ^{L1} )  \text{ has an Inverse-G-Wishart} (G_{full}, \xi_{q( \Sigma ^{L1})}, \Lambda_ {q( \Sigma ^{L1})} ) \text{ density function} \\
q^{*}(  \Sigma ^{L2} )  \text{ has an Inverse-G-Wishart} (G_{full}, \xi_{q( \Sigma ^{L2})}, \Lambda_ {q( \Sigma ^{L2})}) \text{ density function}
\end{array}
\right.
\label{CH12X_eq_28}
\end{equation*}
and Appendix B.2 of \citeA{nolanwand2020}  provides expressions for the $q$-densities for na\"{i}ve \index{variational inference!mean field variational Bayes approximation (MFVB)}mean field variational Bayesian inference for the parameters (see also \citeA{Lee2016}). But, due to the nature of the random effects design matrices,  $Z^{L1}$ and $Z^{L2}$ are large sparse matrices. Indeed, $X$, $Z^{L1}$ and $Z^{L2}$ are respectively $(N \times p)$, $(N \times N_{g_1})$ and $(N \times N_{g_2})$ matrices with $N= \sum_{i=1}^{m} \sum_{j=1}^{n_i} o_{ij}$, $N_{g_1}= \sum_{i=1}^{m} \sum_{j=1}^{n_i} q_1$ and $N_{g_2}= \sum_{i=1}^{m} \sum_{j=1}^{n_i} q_2$. The total number of elements in $[ X, Z^{L1}, Z^{L2} ]$ is $N(p+N_{g_1}+N_{g_2})$. The number of non-zero elements in $[ X, Z^{L1}, Z^{L2} ]$ is $N(p +q_1 + q_2)$. And so by difference, the number of zeros elements in $[ X, Z^{L1}, Z^{L2} ]$ is $N(N_{g_1}+N_{g_2} - q_1 - q_2)$. If we go back to our example of $m=20$ {\it arrondissements} and ($n_1 = \cdots  = n_m = 4$) {\it quartiers}, assuming a model with only one explanatory variable and one intercept ($p=q_1=q_2=2$) and assuming $100$ transactions per  {\it quartier} ($o_{ij} = 100, \forall i,j $), then the  total number of elements in $[ X, Z^{L1}, Z^{L2} ]$ is $2,576 \times 10^3$ and the ratio of zero elements is $98.136\%$. \citeA{nolanwand2020} and \citeA{nolan2020} have proposed streamlined variational Bayes approximations with specific solutions to \index{variational inference!three-level sparse matrix least squares problem}multilevel sparse matrix problems using only the non-zero elements of $[ X, Z^{L1}, Z^{L2} ]$ (in our case using only $1.864\%$ of all the elements). They show that the following relatively small sub-blocks of $\Sigma_{q(\beta , u)}$ are required for variational inference concerning $\sigma^{2}_{\varepsilon}$, $\Sigma ^{L1}$ and $\Sigma ^{L2}$:
\begin{equation}
\begin{array}{l}
\Sigma_{q(\beta)}, \Sigma_{q(u^{L1}_{i})}, \Sigma_{q(u^{L2}_{ij})}, E_{q} \left[ \left( \beta - \mu_{q(\beta)}  \right)   \left(  u^{L1}_{i} -  \mu_{q( u^{L1}_{i})}     \right)^{\prime}  \right] \\
E_{q} \left[ \left( \beta - \mu_{q(\beta)}  \right)   \left(  u^{L2}_{ij} -  \mu_{q( u^{L2}_{ij})}     \right)^{\prime}  \right] \text{ and } E_{q} \left[ \left(    u^{L1}_{i} -  \mu_{q( u^{L1}_{i})}          \right)   \left(  u^{L2}_{ij} -  \mu_{q( u^{L2}_{ij})}     \right)^{\prime}  \right] 
\end{array}
\label{CH12X_eq_29} 
\end{equation}
and the \index{variational inference!mean field variational Bayes approximation (MFVB)}MFVB updates of $\mu_{q(\beta)}$ and each of the sub-blocks of $\Sigma_{q(\beta , u)}$ in (\ref{CH12X_eq_29}) are expressible as a \index{variational inference!three-level sparse matrix least squares problem}three-level sparse matrix least squares problem. \citeA{nolanwand2020} propose a QR-decomposition-based streamlined algorithm for obtaining MFVB approximate posterior density functions for the parameters in the \index{variational inference!three-level linear mixed model}three-level linear mixed model.\footnote{For a balanced model, with $n_i=n$ and $o_{ij}=o$, \citeA{nolanwand2020} show that the order of magnitudes of the number of floating point operations for na\"{i}ve MFVB and QR-decomposition-based streamlined MFVB algorithms are $O(m^3n^3o)$ and $O(mno)$ and therefore the streamlined MFVB can be $78$ times to $2,300$ times faster than the na\"{i}ve MFVB depending on the size of $m$. The R code of this QR-decomposition-based streamlined algorithm is available at http://matt-wand.utsacademics.info/statsPapers.html.}\\
As an illustration, we apply \citeA{nolanwand2020}'s method on a subset of the \index{housing!model example!Paris}housing transaction prices data of \citeA{Baltagi2015}. We use only the year 2003 and the associated transaction prices of the $10,983$ sold flats. As shown in Figure \ref{fig:figure_paris}, Paris is divided into $m=20$ {\it arrondissements}, each of which is divided into $n_i=n=4$ {\it quartiers} (denoted from 1 to 80 on Figure \ref{fig:figure_paris}).\footnote{However, for the writing of the three-level model and for the Figures, we re-number the {\it quartiers}  from 1 to 4. Thus, the {\it quartiers} (5,6,7,8) of the 2nd {\it arrondissement} are noted (1,2,3,4), ... and the {\it quartiers} (77, 78, 79, 80) of the 20th {\it arrondissement} are also noted (1,2,3,4).}

 %------------------------------- Figure -------------------------------------------------
\begin{figure}[ht]
\centering
\includegraphics[scale=0.30]{Figures/Paris_map.png}
\caption{Administrative division of Paris.}
\label{fig:figure_paris}
\end{figure}
%--------------------------------------------------------------------------------
To fit the logarithm of the housing price (in euros per square meter), we use 6 covariates: the surface of the flat (in sq.m), the number of bathrooms, the garage plots, the number of maid's rooms, the floor level of the flat and the presence (or not) of a balcony.
\begin{equation*}
\begin{array}{lll}
y_{ij}  & = & \sum_{k=1}^{p=7} \left( \beta_k + u^{L1}_{i,k} + u^{L2}_{ij, k} \right) X_{ij,k} + \varepsilon_{ij}   \\
& &   1 \leq i \leq m=20  \text{, } 1 \leq j \leq n_i (\equiv n) =4   \text{, } \forall i 
\end{array}
 \label{CH12X_eq_30} 
\end{equation*}
where $\left( \beta_1 + u^{L1}_{i,1} + u^{L2}_{ij, 1} \right)$ is the random intercept. Results for the fixed effects estimates ($\mu_{q(\beta)}$) and their standard errors (from $\Sigma_{q(\beta)}$) are given in Table \ref{tab:estim_FE}. Coefficients are significantly different from zero at $5\%$ level and this model fits well the data (the $R^2$ is close to one). Estimate of $\sigma_{\varepsilon}$ ($=0.259$) is small as compared to the mean $(=8.194)$ of the dependent variable.
But more interestingly, while the estimated variance of the fixed effect of the intercept ($\Sigma_{q(\beta_1)}$) is small, the mean of the estimated variances of the random effects of the outer group ($ (1/m) \sum_{i=1}^{m} \Sigma_{q(u^{L1}_{i,1})} $) is larger while the mean of the estimated variances of the random effects of the inner group is smaller ($ (1/(m.n)) \sum_{i=1}^{m} \sum_{j=1}^{n}\Sigma_{q(u^{L2}_{ij,1})} $). Thus, there is indeed more heterogeneity between {\it arrondissements} than between  {\it quartiers}  within the same {\it arrondissement}. 
\begin{table}[ht!]
\centering
\begin{tabular}{lrrrr}
  \hline
 & $\mu_{q(\beta)}$ & s.e $\left( \Sigma_{q(\beta)} \right)$ & $t\_$value & Pr($>$$|$t$|$) \\ 
  \hline
intercept & 8.02490 & 0.02778 & 288.91765 & 0.0 \\ 
  surface & 0.00122 & 0.00048 & 2.54670 & 0.01089 \\ 
  bathrooms & 0.05808 & 0.02171 & 2.67596 & 0.00746 \\ 
  garage plots & 0.04377 & 0.02358 & 1.85599 & 0.06348 \\ 
  maid rooms & 0.08867 & 0.02807 & 3.15910 & 0.00159 \\ 
  floor levels & 0.01214 & 0.00616 & 1.97215 & 0.04862 \\ 
  balcony & 0.15424 & 0.02996 & 5.14748 & 0.0 \\ 
   \hline
 & statistics &   &    &    \\
   \hline
 N & 10,983 &   &    &    \\
 $\sigma^{2}_{\varepsilon}$  &  0.06704  &   &    &    \\
  $R^2$ & 0.99901  &   &    &    \\
 $\log \underline{p}\left( y,q\right)$  & -1277.57706  &   &    &    \\
  $\Sigma_{q(\beta_1)}$ & 0.00077  &   &    &    \\
 $ (1/m) \sum_{i=1}^{m} \Sigma_{q(u^{L1}_{i,1})} $ & 0.00773  &   &    &    \\
 $ (1/(m.n)) \sum_{i=1}^{m} \sum_{j=1}^{n}\Sigma_{q(u^{L2}_{ij,1})} $ & 0.00127  &   &    &    \\
   \hline
\end{tabular}
\caption{\label{tab:estim_FE}MFVB approximation of the fixed effects parameters of the three-level linear mixed model.}
\end{table}



One advantage of this flexible model is that it provides marginal effects $\left( \beta_k + u^{L1}_{i,k} + u^{L2}_{ij, k} \right)$ for a specific covariable $X_k$, that take into account multi-dimensional random effects. Figure \ref{fig:figure_densities} gives the nonparametric densities of the marginal effect for each of the covariates, both for the outer group $\left( \beta_k + u^{L1}_{i,k}  \right)$ and for the inner group $\left( \beta_k + u^{L1}_{i,k} + u^{L2}_{ij, k} \right)$. The density curves of the marginal effects for the outer group tend to be smoother and more centered around the fixed effect value $\left( \beta_k \right)$ than those of the marginal effects of the inner group.  These curves highlight the differences in heterogeneity of marginal effects at the {\it arrondissement}  and {\it quartier} levels. Thus, an increase of one unit in the number of bathrooms in the flat increases the price per square meter by an average of $5.8\%$ but with a relatively high variability between $3.5\%$ and $7.5\%$. Similarly, the presence or absence of a balcony in the flat increases the price per square meter by an average of $15.4\%$ but with a relatively high variability between $13.5\%$ and $17\%$. But it is especially the 3D histograms of marginal effects that reveal the heterogeneity of such effects across {\it arrondissements}  and {\it quartiers}  (see Figure \ref{fig:figure_hist3D}). These figures clearly show the differences even between  {\it quartiers} in the same {\it arrondissement}. These \index{variational inference!three-level random effects model}three-level random effects models and the associated streamlined MFVB approximation can be useful tools to estimate multi-dimensional housing prices models. In order to gain even more flexibility, we can hope that in the future four-level or even five-level random effects models will be developed in order to induce, for instance, a nested structure as {\it year},  {\it arrondissement}, {\it quartier}, {\it block} and {\it flats} in order to be able to estimate precise multi-dimensional panels.


 %\begin{center}
%\framebox[1.1\width]{Put Figure \ref{fig:figure_densities} here.} 
%\end{center}

% \begin{center}
%\framebox[1.1\width]{Put Figure \ref{fig:figure_hist3D} here.} 
%\end{center}


%------------------------------- Figure -------------------------------------------------
\begin{figure}[htp]
\centering
\begin{tabular}{cc}
\includegraphics[scale=0.150]{Figures/density_surface.png} & \includegraphics[scale=0.150]{Figures/density_bathrooms.png} \\
\includegraphics[scale=0.150]{Figures/density_garage.png} & \includegraphics[scale=0.150]{Figures/density_maid.png} \\
\includegraphics[scale=0.150]{Figures/density_floor.png} & \includegraphics[scale=0.150]{Figures/density_balcony.png}
\end{tabular}
\caption{MFVB approximation of densities of marginal effects. Vertical line: fixed effect $\beta$. Solid line: random effect - {\it arrondissement} $\beta + u^{(L1)}$. Dashed line: random effect - {\it arrondissement + quartier} $\beta + u^{(L1)} + u^{(L2)}$.}
\label{fig:figure_densities}
\end{figure}
%--------------------------------------------------------------------------------


%------------------------------- Figure -------------------------------------------------
%\begin{figure}[htp]
%\centering
%\begin{tabular}{cc}
%\includegraphics[scale=0.150]
%{Figures/dens_surface.png} & %\includegraphics[scale=0.150]
%{Figures/dens_bathrooms.png} \\
%\includegraphics[scale=0.150]
% {Figures/dens_garage.png} & %\includegraphics[scale=0.150]
%{Figures/dens_maid.png} \\
%\includegraphics[scale=0.150]
%{Figures/dens_floor.png} & %\includegraphics[scale=0.150]
%{Figures/dens_balcony.png}
%\end{tabular}
%\caption{MFVB approximation of densities of %marginal effects.}
%\label{fig:figure_densities}
%\end{figure}
%--------------------------------------------------------------------------------
\newpage


%------------------------------- Figure -------------------------------------------------
\begin{figure}[htp]
\centering
\begin{tabular}{ll}
\includegraphics[scale=0.30]{Figures/hist3D_surface.png} & \includegraphics[scale=0.30]{Figures/hist3D_bathrooms.png} \\
\includegraphics[scale=0.30]{Figures/hist3D_garage.png} & \includegraphics[scale=0.30]{Figures/hist3D_maid.png} \\
\includegraphics[scale=0.30]{Figures/hist3D_floor.png} & \includegraphics[scale=0.30]{Figures/hist3D_balcony.png}
\end{tabular}
\caption{3D histograms of marginal effects.}
\label{fig:figure_hist3D}
\end{figure}
%--------------------------------------------------------------------------------


\section{Conclusion}
\label{CH12X_section7}

The development of modelling housing in multi-dimensional frameworks (3D, 4D
or more) is still in its infancy. As compared to the huge literature in a 2D
framework, which explains why there are relatively few multi-dimensional
housing studies. The limitation comes from the availability of the data and
the complexity of the methods relative to time series or longitudinal
dimensions. The previous papers show that both spatial and temporal
dimensions in dynamic systems should be included for \index{hedonic!housing model}hedonic housing models
and discrete models of \index{residential!location}residential location in a \index{three-dimensional}three-dimensional
framework. But the inclusion of these multiple dimensions greatly
complicates the specification and modeling of such systems. Extending models
with unobserved \index{neighborhood}neighborhood characteristics to deal with the endogenous
\index{neighborhood}neighborhood characteristics or introducing rationing in housing markets
(see \citeA{Geyer2013}) is not trivial.

Part of the attractiveness of a \index{neighborhood}neighborhood may be driven by the
characteristics of \index{neighbor}neighbors (for instance, higher-income households attract
higher-income households while lower-income households repel higher-income
households). As \citeA{Kuminoff2013} said \textit{\textquotedblleft
households \textquotedblleft sort" across \index{neighborhood}neighborhoods according to their
wealth and their preferences for public goods, social characteristics, and
commuting opportunities ... These \textquotedblleft equilibrium
\index{sorting}sorting\textquotedblright\ models use the properties of market equilibria,
together with information on household behavior, to infer structural
parameters that characterize preference \index{heterogeneity}heterogeneity. These results can be
used to develop theoretically consistent predictions for the welfare
implications of future policy changes. Analysis is not confined to marginal
effects or a partial equilibrium setting. Nor is it limited to prices and
quantities... These capabilities are just beginning to be understood and
used in applied research"} (p. 1007). Since three decades, econometric
methods have made significant progress and considerably grew to eliminate
non credible assumptions such as \index{homogenous preferences}homogenous preferences and exogenous  \index{amenities}amenities.
But now, in a 2D framework, the structural estimators
still rely on parametric assumptions for utility functions, on specific
statistical distributions (log-normal, Type I extreme value, \index{generalized!extreme value distribution}generalized
extreme value, etc. ) used to capture sources of unobserved \index{heterogeneity}heterogeneity
and some strong assumptions to eliminate potential sources of market
frictions. As suggested by \citeA{Kuminoff2013}, one approach could be to
refine the current estimators through the lens of the econometric literature
on \index{partial identification}partial identification (see \citeA{Manski2007}) which views economic
models as sets of assumptions, some of which are plausible and some of which
are \textquotedblleft esoteric" (according to \citeA{Tamer2010}'s expression)
and are needed only to complete a model. One of the key advantages of such
approach is that it could characterize the potential sensitivity of outcomes
to the least credible assumptions. Last, we have shown that  \index{variational inference!variational Bayesian approximations}variational Bayesian approximations
are  promising future pathways to potentially overcome many problems in applied modelling of housing, hoping that four-level and five-level random effects models will be available soon.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



\section*{Appendix: \index{variational inference!Inverse G-Wishart distributions}Inverse G-Wishart distributions}
\label{CH12X_Appendix}

Inverse G-Wishart distributions are a generalization of inverse Wishart distributions. They correspond to the matrix inverses of random matrices that have a G-Wishart distribution (see \citeA{atay2005monte}, \citeA{nolanwand2020} or \citeA{maestrini2021inverse}). For any positive integer $d$, let $G$ be an undirected graph with $d$ nodes labeled $1, \cdots,d$ and set $E$ consisting of sets of pairs of nodes that are connected by an edge. We say that the symmetric $(d \times d)$ matrix $M$ respects $G$ if $M_{ij}=0$ for all $\left\{i, j \right\} \notin E$. A $(d \times d)$ random matrix $X$ has an Inverse G-Wishart distribution with graph $G$ and parameters $\xi > 0$ and symmetric $(d \times d)$ matrix $\Lambda$, written $X \sim \text{Inverse-G-Wishart}(G,\xi, \lambda)$, if and only if the density function of $X$ satisfies
$$
p \left( X  \right) \propto | X |^{-(\xi+2)/2}\exp \left\{ - \frac{1}{2} \text{tr} (\Lambda X^{-1})   \right\}
$$
 over arguments $X$ such that $X$ is symmetric and positive definite and $X^{-1}$ respects $G$. Two important special cases are $G = G_{full} (\equiv$ totally connected $d$-node graph)
for which the Inverse G-Wishart distribution coincides with the ordinary inverse Wishart distribution, and $G = G_{diag} (\equiv$ totally disconnected $d$-node graph), for which the Inverse G-Wishart distribution coincides with a product of independent inverse
Chi-Squared random variables. The subscripts of $G_{full}$ and $G_{diag}$ reflect the fact that $X^{-1}$ is a full matrix and $X^{-1}$ is a diagonal matrix in each special case.
The $G =G_{full}$ case corresponds to the ordinary inverse Wishart distribution. If $d=1$, the graph $G = G_{full} =G_{diag}$ and the inverse G-Wishart distribution reduces to the Inverse Chi-Squared distribution. 



\bibliography{CH12_Baltagi_Bresson_bibliography}

\printindex


\end{document}