% Author: Simon-Pierre Boucher — contact@spboucher.ai
% =============================================================================
% appendix_methods.tex — Methodological details
% =============================================================================
\section{Methodological Appendix}\label{app:methods}

\subsection{Haversine Distance Formula}\label{app:haversine}

The Haversine formula computes the great-circle distance between two points on the surface of a sphere, given their latitudes and longitudes.  For two points $(lat_1, lon_1)$ and $(lat_2, lon_2)$ expressed in radians, the Haversine distance is:

\begin{equation}\label{eq:haversine}
    d = 2R \cdot \arcsin\!\left(\sqrt{\sin^2\!\left(\frac{\Delta lat}{2}\right) + \cos(lat_1) \cdot \cos(lat_2) \cdot \sin^2\!\left(\frac{\Delta lon}{2}\right)}\right),
\end{equation}

\noindent where:
\begin{itemize}
    \item $R = 6{,}371$ km is the mean radius of the Earth;
    \item $\Delta lat = lat_2 - lat_1$ is the difference in latitudes;
    \item $\Delta lon = lon_2 - lon_1$ is the difference in longitudes.
\end{itemize}

The Haversine formula provides an accurate approximation for short to medium distances on the Earth's surface.  For the distances relevant in our application (typically less than 5 km), the approximation error relative to the exact geodesic distance (Vincenty's formula) is negligible.

\paragraph{Implementation.}  In our spatial merge (Section~\ref{sec:spatial_merge}), the Haversine distance is computed in vectorised form using NumPy.  For each chunk of rental listings, we compute the full pairwise distance matrix between the chunk and all Airbnb listings, yielding an $(n_{\text{chunk}} \times n_{\text{Airbnb}})$ matrix of distances.  Buffer membership is determined by comparing each element of this matrix against the threshold radius $r$.

\subsection{Spatial Weights Construction}\label{app:spatial_weights}

The spatial weights matrix $\mathbf{W}$ used in Model~4 (Equation~\ref{eq:spatial_lag}) is a $k$-nearest-neighbour matrix with $k = 5$.  For each rental listing $i$, the neighbourhood $\mathcal{N}_i$ is the set of the five rental listings closest to $i$ in Euclidean coordinate space (excluding $i$ itself):

\begin{equation}
    w_{ik} = \begin{cases}
        \frac{1}{5} & \text{if } k \in \mathcal{N}_i, \\[4pt]
        0 & \text{otherwise},
    \end{cases}
\end{equation}

\noindent so that the spatial lag $\sum_k w_{ik} \ln(\text{rent}_k)$ is the mean log rent of listing $i$'s five nearest neighbours.  Row-standardisation guarantees that every observation has a well-defined, equally weighted neighbourhood, avoiding the empty-neighbourhood problem that a fixed-radius definition would create in sparse rural areas.

\paragraph{Choice of $k$.}  The choice of $k = 5$ balances two considerations.  A smaller $k$ would capture only the most proximate listings, yielding a responsive but noisy spatial lag; a larger $k$ would average over a broader area, reducing noise but blurring relevant spatial variation---and, in dense urban cores, five nearest neighbours typically lie within a few hundred metres.

\paragraph{Estimation.}  The inclusion of a spatial lag of the dependent variable as a regressor introduces a well-known simultaneity problem: $\sum_k w_{ik} \ln(\text{rent}_k)$ is correlated with $\eta_i$ whenever the errors are spatially correlated, so OLS estimation of Equation~\ref{eq:spatial_lag} would be inconsistent.  We therefore estimate the SAR model by the generalised-moments instrumental-variable estimator of Kelejian and Prucha (\texttt{GM\_Lag} in \texttt{spreg}), which instruments the spatial lag with spatially lagged exogenous regressors, and the SEM by the corresponding \texttt{GM\_Error} estimator \citep{anselin1988spatial, lesage2009introduction}.

\subsection{Quantile Regression Estimation}

The quantile regression model at quantile $\tau$ (Equation~\ref{eq:quantile}) is estimated by minimising the asymmetrically weighted sum of absolute residuals:

\begin{equation}
    \hat{\boldsymbol{\beta}}_\tau = \arg\min_{\boldsymbol{\beta}} \sum_{i=1}^{n} \rho_\tau\!\left(\ln(\text{rent}_i) - \mathbf{X}_i' \boldsymbol{\beta}\right),
\end{equation}

\noindent where $\rho_\tau(u) = u(\tau - \mathbf{1}[u < 0])$ is the check function \citep{koenker1978regression}.  We estimate the model with the iteratively reweighted least squares algorithm implemented in \texttt{statsmodels}.

Standard errors are the asymptotic estimates based on the kernel estimate of the conditional density of the response at the fitted quantile (the default in \texttt{statsmodels}), which are consistent under independent but not necessarily identically distributed errors.

\subsection{Machine-Learning Model Details}

\paragraph{LASSO and Elastic Net.}  The LASSO model \citep{tibshirani1996regression} estimates:
\begin{equation}
    \hat{\boldsymbol{\beta}}_{\text{LASSO}} = \arg\min_{\boldsymbol{\beta}} \left\{ \frac{1}{2n}\|\mathbf{y} - \mathbf{X}\boldsymbol{\beta}\|_2^2 + \lambda \|\boldsymbol{\beta}\|_1 \right\},
\end{equation}
where $\lambda > 0$ is the regularisation parameter selected by 5-fold cross-validation.  The elastic net extends this with an $L_2$ penalty, controlled by a mixing parameter $\alpha \in (0,1)$.

\paragraph{Random Forest.}  The random forest model \citep{breiman2001random} constructs an ensemble of $B = 500$ regression trees with a maximum depth of 15, each trained on a bootstrap sample with random feature subsampling.

\paragraph{Gradient Boosting.}  The gradient boosting model \citep{friedman2001greedy} sequentially fits shallow regression trees to the residuals of the current ensemble.  We use 500 boosting iterations, a maximum tree depth of 5, and a learning rate of 0.05.  SHAP (SHapley Additive exPlanations) values \citep{lundberg2017unified} are computed on the held-out test set from the fitted gradient boosting model to quantify the contribution of each feature to individual predictions.