Geostatistical Analysis of Airbnb Prices

This page is a technical log corresponding to my Airbnb geostatistics project. I made this so I could write about the statistics involved a litle more and explain the code behind the whole approach!

Each code block below is self-contained. Click Run in any section without needing to run the ones before it.

---

Step 1: Simulating Spatially Correlated Data

For the purposes of this simulation, we start with a draw of prices from a multivariate normal distribution where correlation between listings follows an exponential covariance function (Chilès & Delfiner, 2012; Cressie, 1993):

$$Z \sim \mathcal{N}(\mu \mathbf{1},\; \Sigma), \qquad \Sigma_{ij} = \sigma^2 \exp\left(-\frac{|s_i - s_j|}{\phi}\right) + \tau^2 \mathbf{1}_{i=j}$$

Three parameters control the spatial pattern:

• Range $\phi$: How far the correlation extends. Past this, prices are essentially unrelated.
• Partial sill $\sigma^2$: Characterizes how much prices vary due to location.
• Nugget $\tau^2$: Introduces random noise between neighbors.

1a. Set up coordinates

We scatter 200 listings across downtown Austin. The coordinate range (-97.8 to -97.6 longitude, 30.2 to 30.4 latitude) covers roughly 11 by 14 miles. In a real dataset these would be actual lat/lon coordinates scraped from Airbnb listings, but for our purposes we just use uniform random points to demonstrate the spatial methods.

1b. Build the covariance matrix

The covariance matrix $\Sigma$ encodes every pairwise relationship. Entry $\Sigma_{ij}$ tells us how much listing $i$‘s price tells us about listing $j$‘s price. The diagonal entries ( where $i = j$) are the total variance, which is set to $\sigma^2 + \tau^2 = 500$. Non-diagonal entries decay exponentially with distance: two listings on the same block (0.01 degrees apart) have covariance around $400 \cdot \exp(-0.01 / 0.04) \approx 400 \cdot 0.78 = 312$ while two listings on opposite sides of the city (0.15 degrees apart) have covariance near zero.

With $\phi = 0.04$, correlation drops by 95% at roughly $3\phi = 0.12$ degrees, or about 8 miles. At Austin’s latitude (30°N), one degree of longitude is about 0.87 times the distance of one degree of latitude.

Note: Using raw Euclidean degrees distorts the covariance structure at city scale, but this is ignored for this demonstration. This problem is usually handled by using projected coordinates to keep true distances.

1c. Draw prices

Now we sample from the multivariate normal. Because of the spatial covariance structure, nearby listings will tend to have similar prices. Listings far apart can differ by $50 or more even though they share the same mean. This is the spatial structure we’re trying to detect and exploit.

---

Step 2: Sample Variogram

The variogram measures dissimilarity over distance (Cressie, 1993; Matheron, 1963). Take all pairs of points roughly $h$ apart; if their prices are similar, the variogram is low; if prices vary wildly, it’s high:

$$\gamma(h) = \frac{1}{2} \mathbb{E}\left[(Z(s) - Z(s+h))^2\right]$$

The plot should rise from a small value (nearby points are similar) toward a flat plateau (where far apart = uncorrelated).

2a. Bin pairs by distance

The variogram replaces the smooth covariance function with an empirical estimate. Take all $\binom{200}{2}$ listing pairs, group them by how far apart they are, and average the squared price difference within each bin. If spatial structure exists, nearby pairs should have small differences (low semivariance) and distant pairs should have large differences (high semivariance), eventually leveling off at the total variance.

The bin width matters. Too narrow and each bin has few pairs so the data is noisy and prone to variation. Too wide and you can inflate the nugget estimate. Here we use 0.008 degrees (about 0.3 miles), which gives about 15 bins across the 0.15-degree range.

2b. Plot the variogram

Dollar differences between listing pairs organized by distance:

The variogram rises from ~100 (the nugget — noise between neighbors) and approaches ~500 (the sill — total variance). The distance at which it levels off corresponds to the range parameter $\phi = 0.04$.

---

Step 3: Fitting a Variogram Model

The sample variogram is noisy — we need a smooth function for prediction. The exponential model captures the shape with three parameters (Chilès & Delfiner, 2012; Cressie, 1993):

$$\gamma(h) = c_0 + c_1\left(1 - \exp\left(-\frac{h}{a}\right)\right)$$

• Nugget $c_0$ — noise between neighbors
• Partial sill $c_1$ — spatial variance
• Range $a$ — how fast correlation decays. At $3a$, the variogram reaches 95% of the sill.

3a. Fit with scipy

The true parameters we used in Step 1 were nugget = 100, partial sill = 400 (sill = 500), range = 0.04. The fitted values should be close. The nugget is the hardest to properly estimate because it uses only the closest pairs, which are the least common.

3b. Plot sample vs fitted curve

The fitted curve passes through the middle of the noisy sample points. The resulting parameters (nugget ≈ 100, sill ≈ 500, range ≈ 0.04) match the true values we used in Step 1, so our model works!

---

Step 4: Ordinary Kriging

Now we use our model for further price prediction.

Kriging takes a weighted average of observed prices, with weights from the variogram model (Chilès & Delfiner, 2012; Cressie, 1993). Nearby points get more weight, but clustered points share weight so no single group dominates:

$$\hat{Z}(s_0) = \sum_{i=1}^{n} w_i Z(s_i), \qquad \sum_{i=1}^{n} w_i = 1$$

The weights solve a linear system built from the variogram and a Lagrange multiplier $\lambda$ to enforce a sum-to-one constraint:

$$\begin{bmatrix} \Gamma & \mathbf{1} \\ \mathbf{1}^\top & 0 \end{bmatrix} \begin{bmatrix} \mathbf{w} \\ \lambda \end{bmatrix} = \begin{bmatrix} \gamma_0 \\ 1 \end{bmatrix}$$

4a. Build the kriging system

We use 140 training points to predict 60 test points. The true variogram parameters are unknown in actual kriging, but we use them here to focus on how kriging works:

The kriging matrix has a block structure.

The top-left block contains the variogram matrix, $\Gamma = [\gamma(p_i, p_j)]$, where each entry is the variogram value between a pair of training locations. These values encode the spatial relationships within the known data.

The last row and column are ones, $\mathbf{1}$ (with a zero in the bottom-right corner), enforcing the Lagrange constraint we mentioned earlier that the kriging weights satisfy $\sum_{i=1}^{n} w_i = 1$.

This is a general structure. Regardless of the number of training points, the matrix is always of size $(n+1)\times(n+1)$ and is symmetric, with this same block form,

$$\mathbf{K}= \begin{bmatrix} \Gamma & \mathbf{1}\\ \mathbf{1}^\top & 0 \end{bmatrix}.$$

4b. Predict at each test location

For each test point, the right-hand side vector $\gamma_0$ measures its variogram distance to every training point. Solving the kriging system produces one weight for each training point. Each weight is then multiplied by the corresponding known training price, and the weighted prices are summed to obtain the predicted value.

$$\hat{y} = \sum_{i=1}^{n} w_i y_i.$$

An RMSE of about 20 means the typical prediction error is $20. The simulation's total variance is 500 (a standard deviation of about$22), so an RMSE of $20 corresponds to an$R^2$of roughly$1 - 20^2/500 = 0.2$. Thus, we can say using ONLY location captures about a fifth of the variation in Airbnb prices at this scale. Pretty good!

4c. Plot predictions vs actual

Points cluster near the 1:1 line, meaning kriging captures the spatial structure. The scatter around the line reflects the nugget noise. Even if we have perfect spatial knowledge, the individual listings will vary randomly.

---

This analysis is based on material from Nicolas Christou’s UCLA course Statistics C173/C273: Applied Geostatistics. See the project page for the full R-based analysis on real Airbnb data.