% ==========================================================================
\section{Methodology}
\label{sec:methodology}
% ==========================================================================

\subsection{Hedonic Pricing Framework}

Following \citet{Rosen1974}, we model the nightly listing price as a
function of the listing's observable characteristics. Let $p_i$ denote the
nightly price of listing $i$ and let $\mathbf{x}_i$ denote the vector of
$K$ characteristics. The hedonic price function is:
\begin{equation}
p_i = P(\mathbf{x}_i) + \eps_i,
\label{eq:hedonic_general}
\end{equation}
where $P(\cdot)$ maps the characteristic space into prices and
$\eps_i$ is an idiosyncratic error term.

We adopt a log-linear specification:
\begin{equation}
\ln p_i = \alpha + \boldsymbol{\beta}' \mathbf{x}_i + \eps_i,
\label{eq:loglinear}
\end{equation}
where $\alpha$ is a constant, $\boldsymbol{\beta}$ is the $K \times 1$
vector of hedonic coefficients, and $\eps_i \sim (0, \sigma^2_i)$ is
mean-zero but potentially heteroskedastic. The log transformation
accommodates the right-skewed price distribution, yields coefficients
interpretable as approximate percentage changes, and imposes a
multiplicative structure consistent with proportional interactions among
characteristics.

\subsection{Model Specification}

The full hedonic equation expands the characteristic vector into seven
attribute blocks:
\begin{equation}
\begin{aligned}
\ln p_i = \; & \alpha
  + \boldsymbol{\beta}_S' \, \mathbf{S}_i
  + \boldsymbol{\beta}_A' \, \mathbf{A}_i
  + \boldsymbol{\beta}_H' \, \mathbf{H}_i
  + \boldsymbol{\beta}_R' \, \mathbf{R}_i \\
  & + \boldsymbol{\beta}_L' \, \mathbf{L}_i
  + \boldsymbol{\beta}_N' \, \mathbf{N}_i
  + \boldsymbol{\beta}_C' \, \mathbf{C}_i
  + \eps_i,
\end{aligned}
\label{eq:full_model}
\end{equation}
where $\mathbf{S}_i$ collects structural characteristics,
$\mathbf{A}_i$ contains amenity indices,
$\mathbf{H}_i$ captures host attributes,
$\mathbf{R}_i$ comprises reputation variables,
$\mathbf{L}_i$ includes locational controls,
$\mathbf{N}_i$ represents text-derived features, and
$\mathbf{C}_i$ denotes spatial context controls.

\subsection{Feature Engineering}

\textit{Structural features.} We include linear terms for bedrooms,
bathrooms, and guest capacity, plus log-transformed counterparts
($\log(1 + x)$). Per-guest ratios address the proportionality of space
to capacity. Quadratic terms (accommodates$^2$, bedrooms$^2$,
bathrooms$^2$) and interactions (accommodates $\times$ bedrooms,
accommodates $\times$ bathrooms) capture non-linear and joint effects.

\textit{Amenity features.} Six expert-weighted composite indices are
constructed from 24 tracked amenity categories, with weights reflecting
assumed market value (e.g., pool receives weight 2.0 in the luxury index;
Wi-Fi receives weight 3.0 in the remote work index). All indices are
normalized to $[0, 1]$.

\textit{Host features.} Superhost status enters as a binary indicator. A
multi-listing host indicator proxies for professionalization. An interaction
term (Superhost $\times$ entire-home) tests whether the Superhost premium
differs by listing type.

\textit{Reputation features.} The overall review score (0--5) enters
linearly. The number of reviews enters in log form. A composite reputation
score combines normalized rating and review volume:
\begin{equation}
\text{reputation}_i = 0.6 \times \frac{\text{rating}_i}{5}
  + 0.4 \times \frac{\log(1 + \text{reviews}_i)}
                     {\max_j \log(1 + \text{reviews}_j)}.
\label{eq:reputation}
\end{equation}

\textit{Locational features.} Distance to the nearest city centre is
computed using the Haversine formula:
\begin{equation}
d_{ij} = 2R \arcsin\!\left(
  \sqrt{\sin^2\!\left(\frac{\Delta\phi}{2}\right)
  + \cos\phi_i \cos\phi_j \sin^2\!\left(\frac{\Delta\lambda}{2}\right)}
\right),
\label{eq:haversine}
\end{equation}
where $R = 6{,}371$~km. Reference coordinates are provided for 22
Quebec municipalities.

\textit{Spatial context features.} Using a BallTree with haversine metric,
we compute local listing density within 500m and 1~km, local average price
within 500m and 1~km, and a $K$-means spatial cluster ($K = 20$).

\textit{Text features.} From listing names and descriptions, we extract
title length, seven bilingual keyword scores, and 10 TF-IDF/SVD components
(sublinear term frequency, up to 5,000 features, unigrams and bigrams,
dimensionality reduced via Truncated SVD capturing 16.98\% of variance).

\subsection{OLS Estimation with Robust Inference}

We estimate Equation~\eqref{eq:full_model} by OLS with
heteroskedasticity-robust standard errors using the HC3 estimator:
\begin{equation}
\widehat{\Var}(\hat{\boldsymbol{\beta}})_{\text{HC3}}
= (\mathbf{X}'\mathbf{X})^{-1}
  \left( \sum_{i=1}^{N} \frac{\hat{u}_i^2}{(1 - h_{ii})^2}
         \mathbf{x}_i \mathbf{x}_i' \right)
  (\mathbf{X}'\mathbf{X})^{-1},
\label{eq:hc3}
\end{equation}
where $h_{ii}$ is the $i$th diagonal element of the hat matrix
$\mathbf{H} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'$.

\subsection{Nested Model Strategy}

To assess the marginal explanatory power of each attribute block, we
estimate nine nested models:
%
\begin{enumerate}[leftmargin=2em]
    \item[M1.] Structural characteristics only.
    \item[M2.] M1 + amenity bundle.
    \item[M3.] M2 + host attributes.
    \item[M4.] M3 + reputation and reviews.
    \item[M5.] M4 + locational controls.
    \item[M6.] M5 + neighbourhood (spatial cluster).
    \item[M7.] M6 + text-derived features.
    \item[M8.] M7 + spatial context controls.
    \item[M9.] Full hedonic model (identical to M8).
\end{enumerate}
%
We compare models using $R^2$, adjusted $R^2$, AIC, BIC, RMSE, and MAE.

\subsection{Quantile Regression}

To explore heterogeneity across the price distribution, we estimate
quantile regressions \citep{KoenkerBassett1978} at
$\tau \in \{0.10, 0.25, 0.50, 0.75, 0.90\}$:
\begin{equation}
Q_\tau(\ln p_i \mid \mathbf{x}_i)
= \alpha_\tau + \boldsymbol{\beta}_\tau' \, \mathbf{x}_i,
\label{eq:quantile}
\end{equation}
solved by minimizing the asymmetric check function:
\begin{equation}
\min_{\alpha_\tau, \boldsymbol{\beta}_\tau}
\sum_{i=1}^{N} \rho_\tau\!\left(\ln p_i - \alpha_\tau
  - \boldsymbol{\beta}_\tau' \mathbf{x}_i \right),
\quad
\rho_\tau(u) = u \left(\tau - \mathbf{1}(u < 0)\right).
\label{eq:check_function}
\end{equation}
Model fit is assessed via the Koenker--Bassett pseudo-$R^2$:
\begin{equation}
\text{pseudo-}R^2_\tau = 1
  - \frac{\sum_i \rho_\tau(\hat{u}_i)}
         {\sum_i \rho_\tau(\ln p_i - \hat{Q}_\tau^{(0)})}.
\label{eq:pseudoR2}
\end{equation}

\subsection{Machine Learning Models}

We complement the econometric analysis with seven supervised learning
algorithms trained on the same feature set using an 80/20 train/test
split ($N_{\text{train}} = 1{,}524$, $N_{\text{test}} = 382$):
ElasticNet \citep{ZouHastie2005}, Random Forest ($B = 500$ trees)
\citep{Breiman2001}, Gradient Boosting (500 trees, learning rate 0.05)
\citep{Friedman2001}, XGBoost (500 estimators, depth 6)
\citep{ChenGuestrin2016}, LightGBM (500 estimators, 31 leaves),
CatBoost (500 iterations, depth 6), and a Generalized Additive Model.
Predictive performance is evaluated using $R^2$, RMSE, MAE, and MAPE.

\subsection{SHAP Interpretability}

To reconcile predictive power with interpretability, we employ
SHAP \citep{Lundberg2017}. For each prediction $\hat{f}(\mathbf{x}_i)$,
SHAP decomposes the deviation from the baseline into additive feature
contributions:
\begin{equation}
\hat{f}(\mathbf{x}_i) = \phi_0
  + \sum_{k=1}^{K} \phi_k(\mathbf{x}_i),
\label{eq:shap}
\end{equation}
where $\phi_0 = \E[\hat{f}]$ and $\phi_k(\mathbf{x}_i)$ is the SHAP
value of feature $k$ for observation $i$. The mean absolute SHAP value
$\overline{|\phi_k|}$ provides a global importance ranking.

\subsection{Spatial Econometric Models}

We construct a row-standardized $k$-nearest-neighbours weight matrix
$\mathbf{W}$ with $k = 5$ and estimate three spatial models
\citep{LeSagePace2009, Anselin1988}:

\textit{Spatial Lag Model (SLM):}
\begin{equation}
\ln \mathbf{p} = \rho \, \mathbf{W} \ln \mathbf{p}
  + \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\eps},
\label{eq:slm}
\end{equation}
where $\rho$ captures endogenous spatial interaction.

\textit{Spatial Error Model (SEM):}
\begin{equation}
\ln \mathbf{p} = \mathbf{X} \boldsymbol{\beta} + \mathbf{u}, \quad
\mathbf{u} = \lambda \, \mathbf{W} \mathbf{u}
  + \boldsymbol{\eps},
\label{eq:sem}
\end{equation}
where $\lambda$ captures spatial autocorrelation in the errors.

\textit{Spatial Durbin Model (SDM):}
\begin{equation}
\ln \mathbf{p} = \rho \, \mathbf{W} \ln \mathbf{p}
  + \mathbf{X} \boldsymbol{\beta}
  + \mathbf{W} \mathbf{X} \boldsymbol{\theta}
  + \boldsymbol{\eps},
\label{eq:sdm}
\end{equation}
which nests both SLM and SEM under parameter restrictions and allows
neighbours' characteristics to exert direct spillover effects.

Spatial dependence is tested using Moran's $I$ on the OLS residuals:
\begin{equation}
I = \frac{N}{\sum_{i}\sum_{j} w_{ij}}
    \frac{\sum_{i}\sum_{j} w_{ij}
          (\hat{u}_i - \bar{\hat{u}})(\hat{u}_j - \bar{\hat{u}})}
         {\sum_{i} (\hat{u}_i - \bar{\hat{u}})^2}.
\label{eq:moranI}
\end{equation}

\subsection{Price Decomposition}

For listing $i$, the predicted log-price decomposes as:
\begin{equation}
\widehat{\ln p_i}
= \hat{\alpha}
  + \underbrace{\hat{\boldsymbol{\beta}}_S' \mathbf{S}_i}_{\text{property}}
  + \underbrace{\hat{\boldsymbol{\beta}}_A' \mathbf{A}_i}_{\text{amenity}}
  + \underbrace{\hat{\boldsymbol{\beta}}_H' \mathbf{H}_i}_{\text{host}}
  + \underbrace{\hat{\boldsymbol{\beta}}_R' \mathbf{R}_i}_{\text{review}}
  + \underbrace{\hat{\boldsymbol{\beta}}_L' \mathbf{L}_i}_{\text{location}}
  + \underbrace{\hat{\boldsymbol{\beta}}_N' \mathbf{N}_i}_{\text{text}}
  + \underbrace{\hat{\boldsymbol{\beta}}_C' \mathbf{C}_i}_{\text{spatial}}.
\label{eq:decomposition}
\end{equation}
This decomposition provides hosts with a transparent breakdown of how the
market values each attribute of their listing.