Hur påverkar bergvärme priset på en bostad?

Vi undersöker hur bergvärme påverkar värdet på ett hus genom att se hur mycket felet i Boolis statistiska värdering skiljer sig mellan hus med bergvärme och hus utan bergvärme. För mer detaljer kring metoden kan du läsa ett tidigare blogginlägg i ämnet på https://jobb.booli.se/blog/posts/18473-statistiska-varderingar-som-matt-pa-forsaljningsresultat. I korthet går metoden ut på att vi gör statistiska värderingar av en stor mängd bostäder och sedan delar upp bostäderna i olika grupper för att se om grupperna skiljer sig åt.

Våra datakällor består av Boolis databas med information om bostäder och värderingar, samt brunnsarkivet från SGU. Brunnsarkivet finns på https://www.sgu.se/grundvatten/brunnar-och-dricksvatten/brunnsarkivet/ och innehåller information om var i landet det finns grävda brunnar och vad brunnarna används till. Vi har gjort antagandet att om en fastighet har en energibrunn så används den brunnen till bergvärme/-kyla.

Vi börjar med att läsa in datan som består av alla försäljningar med tillhörande statistiska värderingar i Sverige mellan 2016-06-01 och 2019-06-01. Sedan matchar vi de försäljningarna till brunnsdatan.

Inläsning av data

In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Read residence to sold map and evaluation data
sold2residence = pd.read_csv('data/sold-res-all.csv')
evaluations = pd.read_csv('data/knn_all_ot_77104_20160601-20190601.csv')
evaluations = evaluations[evaluations['objectType'].isin(('Villa', 'Radhus', 'Parhus', 'Kedjehus'))]

# Merge datasets to map residenceId to an evaluation
evals_with_residences = evaluations.merge(sold2residence[['id', 'residenceId', 'areaId', 'areaName']],
                                          left_on='booliId', right_on='id', how='inner')

Vi har data om brunnar med fastighetsbeteckning, kopplat till Boolis bostads-ID, uppdelade i en fil per län. Så vi itererar varje fil och slår samman dem.

In [2]:
# Read well data
from pathlib import Path

partial_dfs = []
pathlist = Path('data/brunnar/').glob('*.csv')
for path in pathlist:
    df = pd.read_csv(str(path))
    partial_dfs.append(df)
wells = pd.concat(partial_dfs, axis=0, sort=False)
energy_wells = wells[wells['ANVANDNING'] == 'ENE'].copy()

Brunnsarkivet kan innehålla information om fler än 1 brunn per fastighet. Vi är intresserade av att undersöka om det funnits en energibrunn på fastigheten vid försäljningstillfället, så vi väljer ut den tidigast borrade brunnen och rensar bort information om senare borrade brunnar.

Informationen om vilket datum brunnarna borrades finns inte alltid, och när det finns så skiftar formatet något. Vi får skriva en funktion för att tolka datumen och göra om dem till datetime-objekt.

In [3]:
# Attempt to parse well drill dates
def parse_drill_date(date):
    parsed = np.nan
    if not np.isnan(date):
        datestring = str(int(date))
        if len(datestring) == 4:
            datestring = '{}0101'.format(datestring)
        elif len(datestring) == 6:
            datestring = '{}01'.format(datestring)
        try:
            parsed = pd.to_datetime(datestring)
        except:
            pass
    return parsed

energy_wells['drillDate'] = energy_wells['BORRDATUM'].apply(parse_drill_date)

När vi har datumen på datetime-format kan vi sortera vår dataframe på borrdatum och välja ut den tidigast borrade energibrunnen för varje fastighet.

In [4]:
# Sort by date
energy_wells = energy_wells.sort_values("drillDate", ascending=True)

# Keep only the oldest well for each cadastral
energy_wells = energy_wells.drop_duplicates(subset="cadastral", keep='first')
In [5]:
# Merge datasets
combined = evals_with_residences.merge(energy_wells, on='residenceId', how='left')
In [6]:
display("Antal försäljningar: {}".format(evaluations['booliId'].count()))
display("Antal fastighetsmatchade försäljningar: {}".format(evals_with_residences['booliId'].count()))
display("Antal brunnar: {}".format(wells['BRUNNS_ID'].count()))
display("Antal energibrunnar: {}".format(energy_wells['BRUNNS_ID'].count()))
display("Antal försäljningar efter brunnsmatchning: {}".format(combined['booliId'].count()))
'Antal försäljningar: 247672'
'Antal fastighetsmatchade försäljningar: 227430'
'Antal brunnar: 661628'
'Antal energibrunnar: 286279'
'Antal försäljningar efter brunnsmatchning: 227430'

Vi tittar kort på vår data för att verifiera att allt gått bra. Vi kan se att vårt ursprungliga dataset med försäljningspriser och värderingar var på ca 248 000 observationer av husförsäljningar. Efter en matchning mot fastighetsdata har vi kvar ca 227 000 obs.

I brunnsarkivet finns data om totalt ca 662 000 brunnar, varav ca 286 000 är energibrunnar.

Efter en matchning mellan fastighetsbeteckningarna i brunnsarkivet och Boolis försäljningsdata har vi kvar våra 227 000 obs, men nu kompletterad med information energibrunnar.

Förbehandling och filtrering

När vi har försäljningsdatum och borrdatum kan vi beräkna en binär feature som säger om det fanns en borrad brunn vid försäljningstillfället. Vi kallar kolumnen drilled_before_sale och använder den för vår analys.

In [7]:
# Cast solddate string to datetime object
combined['soldDate'] = pd.to_datetime(combined['soldDate'], format='%Y-%m-%d')

# Assign binary class to objects that had a well when they were sold
combined['has_energy_well'] = (~combined['BRUNNS_ID'].isna()) & (combined['ANVANDNING'] == 'ENE')
combined['drilled_before_sale'] = combined.apply(lambda x: x.soldDate > x.drillDate, axis=1)

Vi gör en grundläggande outlier-filtrering på försäljningspris och värderingsfel.

In [8]:
# Remove outliers based on soldPrice and estimation error
filtered = combined
quantiles = filtered['soldPrice'].quantile([0.01, 0.99])
filtered = filtered[(filtered['soldPrice'] >= quantiles[0.01]) &
                    (filtered['soldPrice'] <= quantiles[0.99])]
quantiles = filtered['error_percent'].quantile([0.05, 0.95])
filtered = filtered[(filtered['error_percent'] >= quantiles[0.05]) &
                    (filtered['error_percent'] <= quantiles[0.95])]
In [9]:
display(filtered.groupby('drilled_before_sale')['drilled_before_sale'].count())
drilled_before_sale
False    182859
True      17769
Name: drilled_before_sale, dtype: int64

Det ger oss totalt ca 200 000 observationer i vårt dataset av husförsäljningar, varav ca 18 000 hade bergvärme när de såldes.

Analys

När all data är på plats analyserar vi hur feldistributionen för hela populationen ser ut. Sedan undersöker vi om feldistributionerna skiljer sig åt mellan de olika grupperna. Vi vill också undersöka om det finns någon korrelation mellan förekomsten av bergvärme och geografi, tid och annan bostadsdata.

In [10]:
filtered['error_percent'].hist(bins=20)
display("Median absolute error: {} %".format(filtered['abs_error_percent'].median() * 100))
display(filtered['error_percent'].describe())
'Median absolute error: 13.566440662406626 %'
count    200628.000000
mean          0.030159
std           0.237572
min          -0.403403
25%          -0.128409
50%          -0.002232
75%           0.145479
max           0.842180
Name: error_percent, dtype: float64

Eftersom vi tittar på samtliga småhus i hela landet får vi, trots vår outlier-filtrering, ganska långa svansar med stora värderingsfel. Upp mot 80 % övervärdering eller 40 % undervärdering i de värsta fallen. Felkurvan är centrerad kring 3 %, och standardavvikelsen är ca 24 %, vilket innebär att vi i genomsnitt övervärderar med 3 % och att 2/3 av alla värderingar hamnar inom 24 % från rätt försäljningspris. Medianen av absolutvärdet på felet är ca 13,6 %.

Felkurvan är inte normalfördelad, utan ser snarare lognormal ut. Vi tittar vidare på hur felkurvorna ser ut för de olika klasserna.

Varierar felfördelningarna mellan de som har bergvärme och de som inte har?

In [11]:
f = plt.figure(figsize=(16,4))
f.subplots_adjust(hspace=0.4, wspace=0.4)

for i, label in enumerate(filtered['drilled_before_sale'].unique()):
    df = filtered[filtered['drilled_before_sale'] == label]
    ax = f.add_subplot(1, 2, i+1)
    ax.set_title("Har bergvärme = {} [SEK]".format(label))
    df['error_percent'].hist(bins=20)
    display("Har bergvärme = {}".format(label))
    print("Median absolute error: {} %".format(df['abs_error_percent'].median() * 100))
    print(df['error_percent'].describe())
'Har bergvärme = False'
Median absolute error: 13.529714285714286 %
count    182859.000000
mean          0.035826
std           0.238823
min          -0.403403
25%          -0.123103
50%           0.002423
75%           0.151988
max           0.842180
Name: error_percent, dtype: float64
'Har bergvärme = True'
Median absolute error: 13.881500000000003 %
count    17769.000000
mean        -0.028152
std          0.215820
min         -0.403359
25%         -0.176134
50%         -0.053261
75%          0.080416
max          0.840634
Name: error_percent, dtype: float64

Felkurvorna för båda klasser påminner mycket om varandra, båda har det lognormala utseendet. Felkurvan för de obs som har bergvärme är dock förskjuten till vänster, med ett medelfel på ca -3 %. Det innebär att bostäderna i den klassen sålts för mer än vår värdering. Medianen av absolutfelet är i samma storlek och standardavvikelserna är liknande i båda klasser.

Hur är bergvärmebrunnarna utspridda över tid?

In [12]:
from pandas.plotting import register_matplotlib_converters

register_matplotlib_converters()

drilled = filtered[filtered['drillDate'] >= pd.to_datetime('1976-01-01')]
drilled['drillDate'].hist(bins=30)
Out[12]:
<matplotlib.axes._subplots.AxesSubplot at 0x14a57eb10>

Borrdatumen för energibrunnar har spikar vid åren kring 2005 och 2017. Men det är en förhållandevis jämn spridning från början av 2000-talet fram till idag.

Finns det en geografisk komponent i var man har bergvärme?

In [13]:
import plotly.express as px
import plotly.graph_objects as go

px.set_mapbox_access_token(open(".mapboxtoken").read())

filtered['has_well'] = filtered['drilled_before_sale'].astype('str')
sthlm = filtered[filtered['areaName'] == 'Stockholms län']
fig = px.scatter_mapbox(sthlm,
                        lat="lat",
                        lon="lng",
                        color="has_well",
                        color_continuous_scale=px.colors.cyclical.IceFire,
                        size_max=0.3,
                        zoom=8)
fig.update_layout(mapbox_style='light')
fig.show()