TeX 52.2%
Python 47.8%
1% ============================================================================2\section{Methodology}3\label{sec:methodology}4% ============================================================================56\subsection{Hedonic Specification}78We adopt the log-linear hedonic model as our baseline:9\begin{equation} \label{eq:hedonic}10\ln P_i = \alpha + \boldsymbol{\beta}' \mathbf{X}_i + \boldsymbol{\gamma}' \mathbf{A}_i + \boldsymbol{\delta}' \mathbf{R}_i + \varepsilon_i11\end{equation}12where $P_i$ is the listing price of property $i$, $\mathbf{X}_i$ is a vector of structural characteristics (bedrooms, bathrooms, stories, log floor area, parking spaces, garage, log land area, log distance to Montreal), $\mathbf{A}_i$ is a vector of amenity dummies (highway, hospital, park, school, transit, golf course proximity), and $\mathbf{R}_i$ is a vector of regional dummies (Montreal, Quebec City, Regions; Gatineau as reference). Continuous variables are centered at their sample means to reduce multicollinearity and facilitate interpretation of the intercept \citep{aiken1991multiple}.1314The log-linear form is chosen on both theoretical and empirical grounds. Theoretically, the multiplicative structure implied by the log transformation reflects the intuition that the value of a housing attribute depends on the overall quality level of the property \citep{halvorsen1981choice}. Empirically, the log transformation reduces right-skewness in the price distribution and yields residuals closer to normality \citep{malpezzi2003hedonic}. We center continuous regressors at their sample means, following \citet{marquardt1980comment}, to reduce the condition number of the design matrix and mitigate numerical multicollinearity without affecting slope coefficient estimates.1516Table~\ref{tab:specification} summarizes the variables included in each segment's specification.1718\input{tables/specification}1920The OLS baseline is estimated with HC3 heteroskedasticity-robust standard errors \citep{long2000using, mackinnon1985some}. The HC3 estimator, which uses leverage-adjusted squared residuals $e_i^2 / (1 - h_{ii})^2$, provides superior finite-sample performance relative to the HC0 and HC1 alternatives, particularly when the design matrix contains high-leverage points \citep{cribari2004asymptotic}. The specification varies by segment: the land segment excludes structural variables (bedrooms, bathrooms, floor area), and the plex segment includes a unit count variable.2122\subsection{Quantile Regression}2324The $\tau$-th conditional quantile of $\ln P_i$ given $\mathbf{X}_i$ is estimated by minimizing the check function:25\begin{equation}26\hat{\boldsymbol{\beta}}(\tau) = \arg\min_{\boldsymbol{\beta}} \sum_{i=1}^{n} \rho_\tau \left( \ln P_i - \mathbf{X}_i' \boldsymbol{\beta} \right)27\end{equation}28where $\rho_\tau(u) = u(\tau - \mathbf{1}(u < 0))$ is the asymmetric absolute loss function \citep{koenker1978regression}. We estimate models at $\tau \in \{0.10, 0.25, 0.50, 0.75, 0.90\}$, spanning the affordable-to-luxury spectrum. The choice of these five quantiles follows the convention in the hedonic literature \citep{zietz2008determinants, liao2012hedonic} and provides sufficient resolution to detect monotonic, U-shaped, or sign-switching patterns in coefficient trajectories.2930The minimization problem is a linear program and can be solved using the simplex method or interior point methods \citep{koenker2005quantile}. We use the Barrodale-Roberts algorithm as implemented in \texttt{statsmodels}, with a maximum iteration count of 5,000 to ensure convergence even for large design matrices. Convergence is monitored and reported for all specifications.3132Model fit is assessed via the pseudo-$R^2$ \citep{koenker1999goodness}:33\begin{equation}34\text{Pseudo-}R^2(\tau) = 1 - \frac{\sum_i \rho_\tau(\hat{u}_i)}{\sum_i \rho_\tau(\ln P_i - \hat{Q}_\tau)}35\end{equation}36where $\hat{Q}_\tau$ is the unconditional sample quantile. This measure compares the weighted sum of absolute deviations from the fitted quantile function to the deviations from the unconditional quantile, analogous to the $R^2$ in mean regression. Values are bounded between 0 and 1, with higher values indicating better fit. Unlike the OLS $R^2$, the pseudo-$R^2$ is not directly comparable across quantiles, as the check function weighting changes with $\tau$ \citep{koenker1999goodness}.3738\subsection{Spatial Weights Matrix}3940We construct a $K$-nearest neighbors ($K = 10$) spatial weights matrix with inverse-distance weighting:41\begin{equation}42w_{ij} = \begin{cases}43d_{ij}^{-1} / \sum_{k \in \mathcal{N}_i} d_{ik}^{-1} & \text{if } j \in \mathcal{N}_i \\440 & \text{otherwise}45\end{cases}46\end{equation}47where $\mathcal{N}_i$ denotes the set of $K$ nearest neighbors of property $i$ based on geographic coordinates, and $d_{ij}$ is the Euclidean distance computed from latitude-longitude pairs. Distances are computed using the Haversine-equivalent flat-Earth approximation, which is accurate for the spatial extent of our study area (approximately 500 km east-west). The KNN approach avoids the island problem that can arise with distance-band weights in markets with irregular spatial density \citep{pace2004spatial, getis2004spatial}.4849The choice of $K = 10$ balances local spatial structure (too few neighbors yield noisy estimates) against global smoothing (too many neighbors dilute local effects). We report sensitivity analyses for $K \in \{5, 15, 20\}$ in Section~\ref{sec:robustness}. The row-normalization of $\mathbf{W}$ ensures that the spatial lag $\mathbf{W}\ln\mathbf{P}$ is a weighted average of neighboring log-prices, facilitating the interpretation of $\rho$ as the elasticity of own price with respect to average neighboring prices.5051The spatial weights matrix is constructed using a KD-tree algorithm for computational efficiency \citep{bentley1975multidimensional}, which provides $O(n \log n)$ nearest-neighbor queries. For our largest segment (land, $n = 9{,}239$), the full $n \times n$ weights matrix is stored in compressed sparse row (CSR) format to minimize memory usage.5253\subsection{Two-Stage Spatial Quantile Regression (2SQR)}5455Following \citet{kim2004quantile}, we augment the quantile regression with an instrumented spatial lag:5657\textbf{First stage (OLS):}58\begin{equation}59\mathbf{W}\ln\mathbf{P} = \mathbf{Z}\boldsymbol{\pi} + \mathbf{v}, \quad \mathbf{Z} = [\mathbf{1}, \mathbf{X}, \mathbf{WX}]60\end{equation}6162\textbf{Second stage (quantile regression):}63\begin{equation} \label{eq:2sqr}64Q_\tau(\ln P_i | \mathbf{X}_i, \widehat{\mathbf{W}\ln P}_i) = \alpha(\tau) + \boldsymbol{\beta}'(\tau) \mathbf{X}_i + \rho(\tau) \widehat{\mathbf{W}\ln P}_i65\end{equation}6667The spatially lagged exogenous variables $\mathbf{WX}$ serve as instruments for the endogenous spatial lag $\mathbf{W}\ln\mathbf{P}$. The exclusion restriction assumes that neighbors' characteristics affect a property's price only through their effect on neighbors' prices, not directly. This assumption is standard in the spatial econometrics literature but may be violated if neighborhood composition effects operate independently of prices \citep{gibbons2015spatial}. The coefficient $\rho(\tau)$ captures the strength of spatial spillovers at quantile $\tau$, and its variation across quantiles reveals whether spatial dependence is itself heterogeneous across the price distribution.6869\subsubsection{Instrument Validity and Weak Instrument Diagnostics}7071The validity of the 2SQR estimator depends on the relevance and exogeneity of the instruments $\mathbf{WX}$. Instrument relevance requires that neighbors' exogenous characteristics are correlated with the average of neighbors' prices---a condition that is almost always satisfied in housing markets, where a neighbor's bedrooms, bathrooms, and floor area directly influence its price. We assess instrument strength using the first-stage $F$-statistic, applying the \citet{stock2002survey} rule of thumb that $F > 10$ indicates adequate instrument strength. We report the first-stage $R^2$ and $F$-statistic for each segment.7273Instrument exogeneity is more difficult to test, as the 2SQR model is exactly identified (the number of excluded instruments equals the number of endogenous regressors in a sense). Following the spatial econometrics convention \citep{kelejian1998generalized}, we assume that the spatial structure fully captures the channel through which neighbors' characteristics affect own prices. This assumption is more plausible when the specification includes a rich set of own-property characteristics and regional fixed effects, which absorb common neighborhood factors.7475\subsubsection{Effect Decomposition}7677In a fully specified SAR model, the marginal effects of covariates are derived from the spatial multiplier matrix $(\mathbf{I} - \rho \mathbf{W})^{-1}$. Following \citet{lesage2009introduction}, the formal decomposition is:78\begin{align}79\frac{\partial \mathbf{y}}{\partial x_{jk}} &= (\mathbf{I} - \rho \mathbf{W})^{-1} \mathbf{I}_n \beta_k80\end{align}81where the average \textit{direct} effect is $n^{-1} \text{tr}[(\mathbf{I} - \rho \mathbf{W})^{-1}] \beta_k$, and the average \textit{indirect} effect is the average row sum of the off-diagonal elements. Computing the full matrix inverse $(\mathbf{I} - \rho \mathbf{W})^{-1}$ is computationally demanding for our sample sizes ($n$ up to $9{,}239$).8283We therefore report \textit{approximate scalar} impacts, which assume a homogeneous spatial structure:84\begin{align}85\widetilde{\text{Total}}_k(\tau) &\approx \frac{\beta_k(\tau)}{1 - \rho(\tau)} \\86\widetilde{\text{Indirect}}_k(\tau) &\approx \frac{\rho(\tau) \cdot \beta_k(\tau)}{1 - \rho(\tau)}87\end{align}8889These scalar approximations are exact when $\mathbf{W}$ is row-normalized and all eigenvalues of $\mathbf{W}$ are equal---a condition that holds approximately for $K$-nearest-neighbor weights with large $n$ \citep{lesage2009introduction}. For our specifications, we verified that the scalar approximation yields values within 3\% of the formal trace-based average direct effects computed on a random subsample. The approximate indirect-to-total ratio equals $\rho(\tau)$, providing an intuitive measure of the contribution of spatial spillovers to total amenity capitalization. Readers should interpret these as approximate magnitudes rather than formally derived LeSage--Pace impacts.9091\subsection{Identification and Interpretation}92\label{sec:identification}9394Several identification and interpretation issues warrant explicit discussion, as they inform the scope of claims that can be supported by our estimates.9596\textbf{Nature of the estimates.} The hedonic model estimates conditional associations in listing prices, not structural demand parameters. In a single cross-section, the hedonic price function reflects a market equilibrium that jointly determines prices and the spatial allocation of households \citep{ekeland2004identification, bajari2005estimating}. The quantile coefficients $\hat{\boldsymbol{\beta}}(\tau)$ are best interpreted as \textit{quantile-specific capitalization gradients}: conditional differences in listing prices associated with amenity proximity at different points of the price distribution.9798\textbf{Amenity endogeneity.} Amenity indicators may be endogenous because households sort into locations based on unobserved preferences that correlate with both amenity access and willingness to pay \citep{bayer2007migration}. For example, households with strong preferences for green space may sort into park-proximate neighborhoods, creating a positive correlation between park proximity and unobserved demand characteristics. Our specification partially addresses this through regional fixed effects and the spatial lag, but we cannot fully rule out sorting-induced endogeneity without quasi-experimental variation in amenity provision.99100\textbf{2SQR identification.} The 2SQR estimator relies on the exclusion restriction that spatially lagged exogenous variables $\mathbf{WX}$ affect the focal property's price only through their effect on neighboring prices. This assumption is standard in the spatial econometrics literature \citep{kelejian1998generalized} but may be violated if neighbors' structural characteristics directly affect the focal property's value---for instance, if a neighbor's large lot generates a visual amenity independent of its price effect. Our rich set of own-property controls and regional fixed effects mitigates but does not eliminate this concern.101102\textbf{Willingness-to-pay interpretation.} The WTP calculations translate semi-log coefficients into dollar-valued implicit prices. These are capitalization measures---the dollar difference in listing price associated with amenity proximity---rather than structural welfare measures or causal willingness-to-pay estimates. The 2SQR estimates should be interpreted as spatially adjusted conditional associations rather than causal effects of amenity access.103104\textbf{Listing prices.} The use of listing prices rather than transaction prices introduces potential measurement error. If the listing-to-sale ratio varies systematically with amenity proximity conditional on the price level, the coefficients may be biased. We discuss this limitation further in Sections~\ref{sec:robustness} and~\ref{sec:discussion}.105106\subsection{Bootstrap Inference}107\label{sec:bootstrap}108109Standard asymptotic inference for quantile regression can be unreliable in the presence of spatial dependence and heteroskedasticity \citep{buchinsky1995estimating}. The asymptotic variance of the quantile regression estimator involves the conditional density of the error term at zero, $f_\varepsilon(0 | \mathbf{x})$, which is difficult to estimate nonparametrically, especially at tail quantiles where data are sparse \citep{koenker2005quantile}.110111We employ the nonparametric paired bootstrap with $B = 200$ replications to compute standard errors for all quantile regression coefficients. Each bootstrap sample draws $n$ observation pairs $(y_i, \mathbf{x}_i)$ with replacement. For each bootstrap replication $b = 1, \ldots, B$:112\begin{enumerate}[nosep]113\item Draw a bootstrap sample of size $n$ with replacement.114\item Estimate the quantile regression at each $\tau$ on the bootstrap sample.115\item Record the coefficient vector $\hat{\boldsymbol{\beta}}^{(b)}(\tau)$.116\end{enumerate}117The bootstrap standard error of $\hat{\beta}_k(\tau)$ is the sample standard deviation across bootstrap replications:118\begin{equation}119\widehat{\text{SE}}_{\text{boot}}[\hat{\beta}_k(\tau)] = \left[ \frac{1}{B-1} \sum_{b=1}^{B} \left( \hat{\beta}_k^{(b)}(\tau) - \bar{\hat{\beta}}_k(\tau) \right)^2 \right]^{1/2}120\end{equation}121122We require at least 50\% of bootstrap replications to converge successfully (i.e., at least 100 out of 200); if fewer converge, we flag the quantile-segment combination as potentially unreliable and fall back to asymptotic standard errors. The choice of $B = 200$ balances computational cost against inference precision. \citet{efron1994introduction} show that $B = 200$ is generally sufficient for standard error estimation (as opposed to confidence interval construction, which may require $B \geq 1{,}000$). We verify that increasing $B$ to 500 does not materially change our conclusions in the robustness section.123124\subsection{Inter-Quantile Heterogeneity Tests}125126To formally test whether implicit prices differ across quantiles, we employ the Wald test:127\begin{equation}128W = \frac{[\hat{\beta}_k(\tau_1) - \hat{\beta}_k(\tau_2)]^2}{\widehat{\text{Var}}[\hat{\beta}_k(\tau_1)] + \widehat{\text{Var}}[\hat{\beta}_k(\tau_2)]} \sim \chi^2(1)129\end{equation}130under the null hypothesis $H_0: \beta_k(\tau_1) = \beta_k(\tau_2)$. The variance in the denominator uses bootstrap standard errors, which are asymptotically valid under regularity conditions \citep{koenker2005quantile}. We note that this test assumes independence between the quantile estimators at $\tau_1$ and $\tau_2$, which is only approximate; a more precise test would account for the covariance between estimators at different quantiles using the joint bootstrap distribution. We test all pairwise combinations $(\tau_1, \tau_2) \in \{(0.10, 0.90), (0.25, 0.75), (0.10, 0.50)\}$ for each amenity variable.131132The omnibus test aggregates individual Wald statistics across all covariates:133\begin{equation}134W_{\text{omnibus}} = \sum_{k=1}^{K} \frac{[\hat{\beta}_k(0.10) - \hat{\beta}_k(0.90)]^2}{\widehat{\text{Var}}[\hat{\beta}_k(0.10)] + \widehat{\text{Var}}[\hat{\beta}_k(0.90)]} \sim \chi^2(K)135\end{equation}136under the joint null that all coefficients are equal between $\tau = 0.10$ and $\tau = 0.90$. This is related to the \citet{koenker1982robust} test for heteroskedasticity, which tests for a specific form of quantile heterogeneity.137138\subsection{Multicollinearity Diagnostics}139140We assess multicollinearity using the variance inflation factor (VIF):141\begin{equation}142\text{VIF}_k = \frac{1}{1 - R_k^2}143\end{equation}144where $R_k^2$ is the $R^2$ from regressing $x_k$ on all other regressors. Following conventional thresholds \citep{obrien2007caution, kutner2004applied}, we consider VIF $> 10$ as indicating problematic multicollinearity and VIF $> 5$ as warranting attention. Centering continuous variables at their sample means eliminates ``non-essential'' multicollinearity (inflated VIFs due to the correlation between linear and polynomial or interaction terms) without affecting the slope estimates \citep{marquardt1980comment}.145146We also report the condition number of the design matrix $\mathbf{X}$, defined as the ratio of the largest to smallest singular value. \citet{belsley1980regression} suggest that condition numbers above 30 indicate potentially harmful multicollinearity.147148\subsection{Spatial Autocorrelation Diagnostics}149150We assess residual spatial autocorrelation using Moran's $I$ statistic:151\begin{equation}152I = \frac{n}{S_0} \frac{\mathbf{z}'\mathbf{W}\mathbf{z}}{\mathbf{z}'\mathbf{z}}153\end{equation}154where $\mathbf{z}$ is the vector of demeaned residuals, $\mathbf{W}$ is the spatial weights matrix, and $S_0 = \sum_{i}\sum_{j} w_{ij}$. Under the null hypothesis of no spatial autocorrelation, the expected value is $E[I] = -1/(n-1)$ and the variance is computed using the analytical formula based on the moments of the weight matrix \citep{cliff1981spatial}. Significance is assessed via the normal approximation $Z = (I - E[I]) / \sqrt{\text{Var}[I]}$. We compute Moran's $I$ for OLS, QR, and 2SQR residuals at each quantile, enabling assessment of whether the spatial lag successfully reduces residual spatial dependence.155