1# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca2"""Rent-based proxies: lease matching, median rents, crude gross rental yields and the user-cost diagnostic.34Data (docs/fundamental_value.md §3). `data/raw/external/rdl.db`, table `rdl_housings`, is a sample of residential5leases (monthly rent, rooms, lease year, optional start date, coordinates) copied from the Immo-Ka platform. Limits6that this module documents but cannot repair: (i) it is a self-selected sample of registered leases, not the rental7stock; (ii) no floor area, so a rooms-to-area bridge to the sales table does NOT exist; (iii) the `rooms` field is an8integer whose meaning the registry does not document (its distribution — mode at 2, 90 % in 1..3 — is consistent with9bedrooms); (iv) `year` is the lease year and `start_date` is missing for about 40 % of rows; (v) the raw table also10contains Ontario cities and implausible rents, which the cleaning below removes.1112Economics. Let R be the annual rent a comparable unit commands and P the price. The one-period user-cost condition13(Poterba 1984; Himmelberg, Mayer & Sinai 2005), in the form of theory.tex eq. (usercost) with the Canadian tax treatment,14 R / P = u = i + τ_p + δ + γ − π^e,15so the fundamental value implied by a rent is V^F = R / u and its ratio to the price is (R/P) / u = gross_yield / u.16Here gross_yield is a *crude* ratio: the median rent of leases with a similar rooms count in the same municipality and17year, annualised, divided by the price of the condo (one unit) or the plex (rent × number of units). It ignores vacancy,18operating costs, quality differences between the rental and owner-occupied stock, and rent regulation; it is a19distribution diagnostic, not a valuation.20"""21from __future__ import annotations2223import json24import sqlite325import urllib.request26from pathlib import Path27from typing import Any2829import numpy as np30import pandas as pd3132from ..config import RAW33from ..evaluation._common import to_pandas3435RDL_DB = RAW / "external" / "rdl.db"36QUARTIER_DB = RAW / "external" / "quartier.db"37BOC_SERIES = "V80691335" # Bank of Canada Valet: "Conventional mortgage, 5-year", weekly, percent (verified 2026-10-27)38BOC_URL = f"https://www.bankofcanada.ca/valet/observations/{BOC_SERIES}/json?start_date=2015-01-01"39BOC_CACHE = RAW / "external" / f"boc_{BOC_SERIES}.json"4041LEASE_RENT_MIN, LEASE_RENT_MAX = 200.0, 15_000.0 # monthly CAD; outside = data-entry error42LEASE_YEAR_MIN, LEASE_YEAR_MAX = 2015, 202643MATCH_MAX_KM = 2.0 # lease → nearest analytic sale (municipality attribution)44ROOMS_BUCKETS = ["0-1", "2", "3", "4+", "all"]45YIELD_CLASSES = ("condo", "plex")46MIN_LEASES_PER_CELL = 104748# user-cost parameters (annual rates). ASSUMPTIONS, documented in docs/fundamental_value.md §4; the mortgage rate is49# the only externally measured input (Bank of Canada series above, annual mean of the weekly observations).50USER_COST_PARAMS: dict[str, dict[str, Any]] = {51 "tau_p": {"value": 0.010, "label": "property tax rate (municipal + school, share of market value)", "source": "assumption; Quebec effective rates cluster around 0.7–1.2 % of market value"},52 "delta": {"value": 0.025, "label": "maintenance and depreciation", "source": "assumption; Himmelberg, Mayer and Sinai (2005) use 2.5 %"},53 "gamma": {"value": 0.020, "label": "risk and liquidity premium", "source": "assumption; Himmelberg, Mayer and Sinai (2005) use 2.0 %"},54 "g_scenarios": {"value": [0.00, 0.02, 0.04], "label": "expected nominal appreciation π^e (three scenarios)", "source": "assumption; bracket of long-run nominal house-price growth"},55}565758# ------------------------------------------------------------------------------------------------ leases5960def rooms_bucket(rooms: pd.Series) -> pd.Series:61 r = pd.to_numeric(rooms, errors="coerce")62 out = pd.Series(pd.NA, index=r.index, dtype="string")63 out[r <= 1] = "0-1"64 out[r == 2] = "2"65 out[r == 3] = "3"66 out[r >= 4] = "4+"67 return out686970def load_leases(path: Path = RDL_DB) -> pd.DataFrame:71 """Read `rdl_housings`, clean, and record the cleaning funnel in `df.attrs['funnel']`.7273 Kept: coordinates present, rent in [200, 15 000] CAD/month, lease year in [2015, 2026], rooms recorded (≥ 0),74 exact duplicates (address, apartment, year, rent, rooms) collapsed to one row.75 """76 con = sqlite3.connect(path)77 df = pd.read_sql_query(78 "SELECT id, full_address, apartment_number, city, zip, lat, lng, price AS rent, rooms, year, start_date, "79 "type_of_accomodation, heating_included, electricity_included, furnishing_included, parking_included, address_slug "80 "FROM rdl_housings", con)81 con.close()82 funnel = {"raw": int(len(df))}83 df = df.dropna(subset=["lat", "lng"])84 funnel["coordinates_present"] = int(len(df))85 df["rent"] = pd.to_numeric(df["rent"], errors="coerce")86 df = df[(df["rent"] >= LEASE_RENT_MIN) & (df["rent"] <= LEASE_RENT_MAX)]87 funnel["rent_in_range"] = int(len(df))88 df["year"] = pd.to_numeric(df["year"], errors="coerce")89 df = df[(df["year"] >= LEASE_YEAR_MIN) & (df["year"] <= LEASE_YEAR_MAX)]90 df["year"] = df["year"].astype(int)91 funnel["year_in_range"] = int(len(df))92 df["rooms"] = pd.to_numeric(df["rooms"], errors="coerce")93 df = df[df["rooms"].notna() & (df["rooms"] >= 0)]94 funnel["rooms_recorded"] = int(len(df))95 df = df.drop_duplicates(subset=["address_slug", "apartment_number", "year", "rent", "rooms"])96 funnel["deduplicated"] = int(len(df))97 df["rooms_bucket"] = rooms_bucket(df["rooms"])98 df["start_date"] = pd.to_datetime(df["start_date"].replace("", pd.NA), errors="coerce")99 df["start_date_present"] = df["start_date"].notna()100 df = df.reset_index(drop=True)101 df.attrs["funnel"] = funnel102 return df103104105def _to_32198(lat: np.ndarray, lng: np.ndarray) -> tuple[np.ndarray, np.ndarray]:106 from pyproj import Transformer107 tr = Transformer.from_crs("EPSG:4326", "EPSG:32198", always_xy=True)108 x, y = tr.transform(np.asarray(lng, dtype=float), np.asarray(lat, dtype=float))109 return np.asarray(x), np.asarray(y)110111112def match_leases(leases: pd.DataFrame, analytic: Any, *, max_km: float = MATCH_MAX_KM, quartier_db: Path | None = QUARTIER_DB) -> pd.DataFrame:113 """Attribute each lease to a municipality (nearest analytic sale within `max_km`, Quebec Lambert EPSG:32198) and,114 when `quartier.db` is available, to a 2021 dissemination area (point in polygon). Leases outside Quebec or farther115 than `max_km` from any sale are left unmatched (`mun_code` NaN). Adds: mun_code, mun_name, region_group,116 match_dist_km, dauid."""117 from scipy.spatial import cKDTree118119 a = to_pandas(analytic)120 cols = ["mun_code", "mun_name", "region_group", "lat", "lng"]121 ref = a[cols].dropna(subset=["lat", "lng"]).drop_duplicates(subset=["lat", "lng"])122 rx, ry = _to_32198(ref["lat"].to_numpy(), ref["lng"].to_numpy())123 tree = cKDTree(np.column_stack([rx, ry]))124 lx, ly = _to_32198(leases["lat"].to_numpy(), leases["lng"].to_numpy())125 dist, idx = tree.query(np.column_stack([lx, ly]), k=1, distance_upper_bound=max_km * 1000.0)126 ok = np.isfinite(dist)127 out = leases.copy()128 out["match_dist_km"] = np.where(ok, dist / 1000.0, np.nan)129 for c in ("mun_code", "mun_name", "region_group"):130 vals = ref[c].to_numpy()131 out[c] = pd.Series(np.where(ok, vals[np.where(ok, idx, 0)], None), index=out.index, dtype="object")132 out["dauid"] = pd.NA133 if quartier_db is not None and Path(quartier_db).exists():134 out.loc[ok, "dauid"] = _point_in_da(out.loc[ok, "lat"].to_numpy(), out.loc[ok, "lng"].to_numpy(), Path(quartier_db))135 out["matched"] = ok136 return out137138139def _point_in_da(lat: np.ndarray, lng: np.ndarray, quartier_db: Path) -> np.ndarray:140 from shapely import STRtree141 from shapely.geometry import MultiPolygon, Point, Polygon142143 con = sqlite3.connect(quartier_db)144 geoms, ids = [], []145 for dauid, poly in con.execute("SELECT dauid, poly FROM da_poly"):146 rings = json.loads(poly)147 try:148 if rings and isinstance(rings[0][0][0], (int, float)):149 g = Polygon(rings[0], rings[1:])150 else:151 g = MultiPolygon([Polygon(p[0], p[1:]) for p in rings])152 except (IndexError, TypeError, ValueError):153 continue154 if not g.is_valid:155 g = g.buffer(0)156 geoms.append(g)157 ids.append(str(dauid))158 con.close()159 tree = STRtree(geoms)160 pts = [Point(float(x), float(y)) for x, y in zip(lng, lat)]161 res = np.full(len(pts), None, dtype=object)162 q = tree.query(pts, predicate="within")163 # several polygons can claim a boundary point: keep the first164 seen: set[int] = set()165 for pi, gi in zip(q[0], q[1]):166 if pi not in seen:167 res[pi] = ids[gi]168 seen.add(int(pi))169 return res170171172def median_rents(leases: pd.DataFrame, *, by: tuple[str, ...] = ("mun_code", "year"), min_n: int = 1) -> pd.DataFrame:173 """Median monthly rent by `by` × rooms bucket, plus an "all" bucket. Columns: *by, rooms_bucket, n_leases,174 rent_median, rent_p25, rent_p75. Only matched leases (mun_code present) are used."""175 df = leases.dropna(subset=["mun_code"]).copy()176 df["mun_code"] = df["mun_code"].astype(str)177 keys = list(by)178 frames = []179 for bucket_df, label in ((df, None), (df.assign(rooms_bucket="all"), "all")):180 g = bucket_df.groupby(keys + ["rooms_bucket"], observed=True)["rent"]181 s = g.agg(n_leases="size", rent_median="median", rent_p25=lambda v: v.quantile(0.25), rent_p75=lambda v: v.quantile(0.75)).reset_index()182 frames.append(s)183 out = pd.concat(frames, ignore_index=True)184 out = out[out["n_leases"] >= min_n]185 out["rooms_bucket"] = out["rooms_bucket"].astype(str)186 return out.sort_values(keys + ["rooms_bucket"]).reset_index(drop=True)187188189# ------------------------------------------------------------------------------------------------ yields190191def gross_yields(analytic: Any, rents: pd.DataFrame, *, classes: tuple[str, ...] = YIELD_CLASSES, condo_bucket: str = "all",192 plex_bucket: str = "all", min_leases: int = MIN_LEASES_PER_CELL) -> pd.DataFrame:193 """Crude gross rental yield for condo and plex sales: 12 × median rent (same municipality, same year, chosen rooms194 bucket) × units / price, where units = 1 for condos and `n_units` for plexes (null → 2). No rooms-to-area bridge195 exists, so the rent is NOT adjusted for the size of the sold unit. Returns one row per matched sale."""196 a = to_pandas(analytic)197 a = a[a["property_class"].isin(list(classes))].copy()198 a["tx_id"] = a["tx_id"].astype(str)199 a["mun_code"] = a["mun_code"].astype(str)200 a["year"] = pd.to_datetime(a["date"]).dt.year.astype(int)201 a["rooms_bucket"] = np.where(a["property_class"] == "condo", condo_bucket, plex_bucket)202 r = rents[rents["n_leases"] >= min_leases][["mun_code", "year", "rooms_bucket", "n_leases", "rent_median"]].copy()203 r["mun_code"] = r["mun_code"].astype(str)204 r["year"] = r["year"].astype(int)205 m = a.merge(r, on=["mun_code", "year", "rooms_bucket"], how="inner", validate="many_to_one")206 units = np.where(m["property_class"] == "plex", pd.to_numeric(m.get("n_units"), errors="coerce").fillna(2).clip(lower=1), 1.0)207 m["units"] = units208 m["annual_rent"] = 12.0 * m["rent_median"] * units209 m["gross_yield"] = m["annual_rent"] / m["price"].astype(float)210 keep = ["tx_id", "date", "year", "mun_code", "mun_name", "region_group", "property_class", "price", "units", "n_leases",211 "rent_median", "annual_rent", "gross_yield"]212 return m[[c for c in keep if c in m.columns]].reset_index(drop=True)213214215# ------------------------------------------------------------------------------------------------ user cost216217def user_cost(i: Any, tau_p: float, delta: float, gamma: float, g: Any) -> Any:218 """u = i + τ_p + δ + γ − π^e (annual rates; `i` and `g` may be arrays)."""219 return np.asarray(i, dtype=float) + tau_p + delta + gamma - np.asarray(g, dtype=float)220221222def implied_value_ratio(gross_yield: Any, u: Any) -> Any:223 """(R / u) / P = gross_yield / u: the rent-implied fundamental value relative to the price."""224 return np.asarray(gross_yield, dtype=float) / np.asarray(u, dtype=float)225226227def fetch_boc_mortgage_rate(cache: Path = BOC_CACHE, *, download: bool = True, timeout: float = 20.0) -> pd.DataFrame:228 """Weekly 5-year conventional mortgage rate (percent) from the Bank of Canada Valet API, cached as JSON.229 Returns an empty frame (never an invented series) when neither the download nor the cache is available."""230 payload: dict[str, Any] | None = None231 if download:232 try:233 with urllib.request.urlopen(BOC_URL, timeout=timeout) as resp: # noqa: S310 (fixed https URL)234 payload = json.loads(resp.read().decode("utf-8"))235 cache.parent.mkdir(parents=True, exist_ok=True)236 cache.write_text(json.dumps(payload), encoding="utf-8")237 except Exception: # network failure → cache238 payload = None239 if payload is None and cache.exists():240 payload = json.loads(cache.read_text(encoding="utf-8"))241 if not payload or "observations" not in payload:242 return pd.DataFrame(columns=["date", "rate_pct"])243 rows = [(o["d"], float(o[BOC_SERIES]["v"])) for o in payload["observations"] if o.get(BOC_SERIES, {}).get("v") not in (None, "")]244 df = pd.DataFrame(rows, columns=["date", "rate_pct"])245 df["date"] = pd.to_datetime(df["date"])246 df.attrs["label"] = payload.get("seriesDetail", {}).get(BOC_SERIES, {}).get("label", "")247 return df.sort_values("date").reset_index(drop=True)248249250def annual_mortgage_rate(weekly: pd.DataFrame, years: list[int]) -> pd.Series:251 """Annual mean of the weekly rate as a decimal; NaN for years without observations (never filled)."""252 if weekly.empty:253 return pd.Series(np.nan, index=years, name="i")254 w = weekly.assign(year=weekly["date"].dt.year)255 m = w.groupby("year")["rate_pct"].mean() / 100.0256 return pd.Series([m.get(y, np.nan) for y in years], index=years, name="i")257