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

\subsection{Hedonic Price Function}
\label{subsec:hedonic}

Let $p_i$ denote the displayed nightly price of listing $i$ and
$\mathbf{x}_i$ the vector of its observed characteristics. Following
\citet{Rosen1974}, the hedonic price function $P(\cdot)$ maps
characteristics into prices, and we approximate it by the semi-logarithmic
form that dominates the housing and short-term rental literatures
\citep{Malpezzi2003, Sirmans2005, WangNicolau2017}:
\begin{equation}
\ln p_i = \alpha
  + \boldsymbol{\beta}_S' \mathbf{S}_i
  + \boldsymbol{\beta}_A' \mathbf{A}_i
  + \boldsymbol{\beta}_H' \mathbf{H}_i
  + \boldsymbol{\beta}_L' \mathbf{L}_i
  + \mu_{c(i)}
  + \eps_i ,
\label{eq:hedonic}
\end{equation}
where $\mathbf{S}_i$ collects the structural variables ($\ln$ bedrooms,
$\ln(1+\text{bathrooms})$, $\ln$ guest capacity, a shared-bathroom
indicator, a studio indicator and five title-category indicators with
\emph{Rental unit} as the reference), $\mathbf{A}_i$ the listing-level
presentation variables ($\ln(1+\text{amenities})$, $\ln(1+\text{photos})$,
pets allowed), $\mathbf{H}_i$ the host and reputation variables (Superhost,
Guest Favourite, $\ln(1+\text{reviews})$, a rated indicator and the
centred rating interacted with it), $\mathbf{L}_i$ the location variables
($\ln(1+\text{km to the nearest centre})$ and $\ln(1+\text{listings within
1 km})$), and $\mu_{c(i)}$ a fixed effect for the city or spatial cell
$c(i)$ in which the listing is located. The error $\eps_i$ is mean zero
conditional on the regressors but is allowed to be heteroskedastic and
spatially correlated within cities.

The semi-logarithmic form has three consequences for interpretation. For
a continuous regressor entered in logarithms the coefficient is an
elasticity. For a continuous regressor in levels, such as the rating,
$100\,\beta$ is the approximate percentage change in price per unit. For
an indicator variable the percentage premium is
\begin{equation}
g = 100\left[\exp\!\left(\hat{\beta} - \tfrac{1}{2}\widehat{\Var}(\hat{\beta})\right) - 1\right],
\label{eq:kennedy}
\end{equation}
the \citet{Kennedy1981} correction of the \citet{HalvorsenPalmquist1980}
transformation, which removes the small-sample bias of
$\exp(\hat{\beta}) - 1$ and is the form we report in the text. The
coefficients in Equation~\eqref{eq:hedonic} are conditional associations
between characteristics and asking prices in a cross-section; they are
not marginal willingness to pay, which would require the stronger
assumptions discussed by \citet{Kuminoff2010} and \citet{Bishop2020}, nor
realised transaction prices.

\subsection{Specifications and Inference}
\label{subsec:specs}

We estimate Equation~\eqref{eq:hedonic} by ordinary least squares in five
nested specifications: M1 contains the structural block only; M2 adds the
presentation block; M3 the host and reputation block; M4 the two location
variables; and M5 the city fixed effects. Cities with at least 15 listings
in the analysis sample receive their own effect; the remaining listings,
which are scattered across several hundred small municipalities, are
grouped into 25 spatial cells obtained by $k$-means clustering of their
projected coordinates, so that every listing belongs to exactly one
fixed-effect group. The fixed effects absorb all between-city variation in
prices, including tourism demand, seasonality of the local market and
regulation, and the remaining coefficients are identified from
within-city variation.

Standard errors are heteroskedasticity-robust (HC3,
\citealp{MacKinnonWhite1985}) in all specifications. Because the
regressors of neighbouring listings are correlated and the error term is
likely to be spatially correlated, we also report standard errors
clustered by fixed-effect group for M5 \citep{CameronMiller2015}, and we
re-estimate M5 with 60 finer spatial cells as the fixed-effect unit in the
robustness section. Multicollinearity among the capacity variables is
monitored through variance inflation factors, and block significance is
assessed with robust $F$ tests.

\subsection{Quantile Regression}
\label{subsec:quantile_methods}

To examine whether implicit prices differ between budget and premium
listings we estimate the conditional quantile function of
\citet{KoenkerBassett1978},
\begin{equation}
Q_\tau(\ln p_i \mid \mathbf{x}_i) = \alpha_\tau + \boldsymbol{\beta}_\tau' \mathbf{x}_i + \mu_{c(i),\tau},
\qquad \tau \in \{0.10, 0.25, 0.50, 0.75, 0.90\},
\label{eq:quantile}
\end{equation}
with the full M5 regressor set, by minimising the check-function loss
$\sum_i \rho_\tau(\ln p_i - \alpha_\tau - \boldsymbol{\beta}_\tau'\mathbf{x}_i - \mu_{c(i),\tau})$,
$\rho_\tau(u) = u\,(\tau - \mathbf{1}\{u<0\})$
\citep{KoenkerHallock2001}. Standard errors use the kernel estimate of the
sparsity function. Goodness of fit is the Koenker--Machado pseudo-$R^2$,
one minus the ratio of the minimised loss to the loss of the
unconditional quantile. Because the coefficients at different quantiles
are estimated from the same sample, differences between the
$\tau = 0.90$ and $\tau = 0.10$ coefficients are tested with a paired
bootstrap of 100 replications in which both quantile regressions are
re-estimated on each resample.

\subsection{Spatial Econometric Models}
\label{subsec:spatial_methods}

Listings cluster in space and their prices may depend on neighbours'
prices or on spatially correlated unobservables. We construct a
row-standardised $k$-nearest-neighbour weights matrix $\mathbf{W}$
($k = 5$, with $k = 10$ as a check) from coordinates projected to
kilometres, so that the anisotropy of degrees of longitude at Quebec's
latitude does not distort distances, and we test the OLS residuals of M4
and M5 for spatial autocorrelation with Moran's $I$,
\begin{equation}
I = \frac{N}{\sum_i \sum_j w_{ij}}\,
    \frac{\sum_i \sum_j w_{ij}\,\hat{u}_i \hat{u}_j}{\sum_i \hat{u}_i^2},
\label{eq:moran}
\end{equation}
with inference from 999 random permutations. We then estimate by maximum
likelihood \citep{Anselin1988, LeSagePace2009} the spatial lag model
(SAR), the spatial error model (SEM) and the spatial Durbin model (SDM),
\begin{align}
\text{SAR:}\quad & \ln \mathbf{p} = \rho\,\mathbf{W}\ln\mathbf{p} + \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\eps}, \label{eq:sar}\\
\text{SEM:}\quad & \ln \mathbf{p} = \mathbf{X}\boldsymbol{\beta} + \mathbf{u}, \qquad \mathbf{u} = \lambda\,\mathbf{W}\mathbf{u} + \boldsymbol{\eps}, \label{eq:sem}\\
\text{SDM:}\quad & \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{align}
using the M4 regressor set, since the city fixed effects and the spatial
lag compete for the same between-city variation. The SDM nests the SAR
($\boldsymbol{\theta} = \mathbf{0}$) and the SEM
($\boldsymbol{\theta} = -\rho\boldsymbol{\beta}$), and we use a
likelihood-ratio test of the first restriction \citep{Elhorst2010}.

In models with a spatial lag the coefficients are not marginal effects.
Writing $\mathbf{S}_k(\mathbf{W}) = (\mathbf{I} - \rho\mathbf{W})^{-1}(\mathbf{I}\beta_k + \mathbf{W}\theta_k)$
for regressor $k$ (with $\theta_k = 0$ in the SAR), the average direct
effect is $N^{-1}\,\mathrm{tr}\,\mathbf{S}_k(\mathbf{W})$, the average
total effect is $N^{-1}\,\boldsymbol{\iota}'\mathbf{S}_k(\mathbf{W})\boldsymbol{\iota}$,
and the indirect (spillover) effect is their difference
\citep{LeSagePace2009}. We compute the three summary measures exactly
from the dense inverse and obtain their standard errors from 300 draws
of the parameters from their estimated asymptotic distribution.

\subsection{Machine Learning Benchmarks and Interpretation}
\label{subsec:ml_methods}

The log-linear form imposes additive, constant implicit prices. To
measure how much predictive accuracy this restriction costs we compare
the M5 specification with an elastic net \citep{ZouHastie2005}, a random
forest \citep{Breiman2001}, gradient boosting \citep{Friedman2001},
XGBoost \citep{ChenGuestrin2016} and LightGBM, all trained on the same
regressors, including the fixed-effect indicators. Performance is
measured out of sample by five-fold cross-validation in two designs:
random folds, which is the usual benchmark, and spatially blocked folds in
which the 20 spatial clusters of listings are assigned to folds as whole
blocks \citep{Roberts2017}. The second design asks whether a model trained
on some areas predicts prices in areas it has never seen, which is the
relevant question for an automated valuation tool and a stricter test of
generalisation when observations are spatially dependent. We report the
out-of-sample $R^2$, root mean squared error and mean absolute error on the
log scale and the median absolute percentage error (MdAPE) on the price
scale, $\operatorname{median}_i |\exp(\hat{y}_i) - p_i| / p_i$. The
first version of this paper reported a mean absolute percentage error
computed on log prices, which is not interpretable as a price error and
is not used here.

The best learner is interpreted with SHAP values \citep{Lundberg2017},
computed exactly for tree ensembles \citep{Lundberg2020}. For a prediction
$\hat f(\mathbf{x}_i)$ the SHAP values $\phi_{ik}$ satisfy
$\hat f(\mathbf{x}_i) = \phi_0 + \sum_k \phi_{ik}$ with
$\phi_0 = \E[\hat f]$, and the mean absolute value
$\overline{|\phi_k|} = N^{-1}\sum_i |\phi_{ik}|$ ranks features by their
average contribution to predicted log prices. SHAP values describe the
fitted model, not the data-generating process; they are attribution
weights for associations captured by the learner and carry no causal
content.

\subsection{Identifying the Nightly-Price Regime}
\label{subsec:mixture_methods}

Section~\ref{sec:data} shows that the displayed price field mixes
nightly prices with total prices of multi-night stays. A price threshold
separates the two regimes imperfectly because the stay total of a cheap
listing can fall below the threshold while the nightly price of a large
chalet can exceed it. We therefore model the regime explicitly with a
two-component mixture of hedonic regressions with common slopes,
\begin{equation}
\ln p_i = \mathbf{x}_i'\boldsymbol{\beta} + \delta\, z_i + \eps_i,
\qquad \eps_i \mid z_i \sim \mathcal{N}(0, \sigma^2_{z_i}),
\qquad \Pr(z_i = 1 \mid \mathbf{w}_i) = \Lambda(\mathbf{w}_i'\boldsymbol{\gamma}),
\label{eq:mixture}
\end{equation}
where $z_i = 1$ indicates the stay regime, $\delta$ is the log of the
typical stay length that the displayed price covers, the regime-specific
variance allows stay lengths to vary across listings, and the regime
probability depends through a logistic link on characteristics
$\mathbf{w}_i$ (Montreal location, review count, capacity, Superhost
status and the hotel and chalet categories) that the descriptive analysis
identifies as predictors of the second mode. The model is estimated by
the EM algorithm, initialised at the \$1,000 rule; the E-step computes
the posterior probability $\hat p_i = \Pr(z_i = 1 \mid \ln p_i,
\mathbf{x}_i, \mathbf{w}_i)$ and the M-step updates
$(\boldsymbol{\beta}, \delta)$ by weighted least squares on the stacked
data, $(\sigma_0, \sigma_1)$ from the weighted residuals and
$\boldsymbol{\gamma}$ by a weighted logistic regression. Listings with
$\hat p_i < 0.5$ constitute the analysis sample used in all subsequent
estimation; the robustness section varies the cut-off, re-weights by
$1 - \hat p_i$ and compares with the simple price-threshold rules.