# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Table of Contents
#
# 1. [Building Causal Models](#Building-Causal-Models)
# 2. [Model Analysis](#Model-Analysis)
#     1. [Causal Effect of Changing Neighborhoods](#Causal-Effect-on-ROI-of-Changing-Neighborhoods---Should-we-choose-North-or-South?)
#     2. [Causal Effect of Different BPBs](#Causal-Effect-on-ROI-of-Different-BPBs)
#     3. [Given we want a certain Bed and Baths, where should we buy?](#Given-we-want-a-certain-Bed-and-Baths,-where-should-we-buy?)
#     4. [Given we want to buy in a certain area, what beds and baths should we buy?](#Given-we-want-to-buy-in-a-certain-area,-what-beds-and-baths-should-we-buy?)
# 3. [Property Examples](#Property-Examples)
# 4. [Next Steps](#Next-Steps)

import pyro
import pyro.distributions as dist
from pyro.infer import Importance, EmpiricalMarginal
from statistics import mean
import torch
import pandas as pd
import numpy as np
import statsmodels.api as sm
from statsmodels.formula.api import ols
import statsmodels.stats.multicomp
pyro.set_rng_seed(101)

df = pd.read_csv("data/cleansed_data.csv")

# # Building Causal Models

# ### Fitting Factor Effects Models for Rent and Zestimate
# - In order to get the parameters for the Rent and Zestimate calculations we use factor effects models(aka ANOVA models)

# Taken from [this](https://pythonfordatascience.org/anova-2-way-n-way/) tutorial

# +
rent_model = ols('Rent ~ C(BPB)', df).fit()

# Seeing if the overall model is significant
print(f"Overall model F({rent_model.df_model: .0f},{rent_model.df_resid: .0f}) = {rent_model.fvalue: .3f}, p = {rent_model.f_pvalue: .4f}")
# -

# Model is signficant

res = sm.stats.anova_lm(rent_model, typ=2)
rent_model_se = (res.sum_sq.Residual / res.df.Residual)**0.5
res

# - Do OLS model for Zestimates
# - Get Neighborhood Marginal
# - Get BPB marginal
# - Get ROI function in the code

# +
zest_model = ols('Zest ~ C(BPB) + C(Neigh)', df).fit()

# Seeing if the overall model is significant
print(f"Overall model F({zest_model.df_model: .0f},{zest_model.df_resid: .0f}) = {zest_model.fvalue: .3f}, p = {zest_model.f_pvalue: .4f}")
# -

res = sm.stats.anova_lm(zest_model, typ=2)
zest_model_se = (res.sum_sq.Residual / res.df.Residual)**0.5
res

# ### Functions for Calculating ROI
# - below are the functions and hyperparameters for calculating ROI

# +
# ROI Consts
INFLATION_RATE = 0.028
MORTGAGE_RATE = 0.036
NUM_YEARS = 15
DOWN_PAYMENT = 0.20


def calculate_monthly(P, mortgage_rate, num_years):
    n = num_years * 12
    monthly_i = mortgage_rate / 12
    numerator = monthly_i * (1 + monthly_i) ** n
    denominator = ((1 + monthly_i) ** n) - 1
    return P * numerator / denominator


def airbnb_income(price, inflation_rate, num_years):
    total = 0
    for year_number in range(num_years):
        curr_inflation = (1 + inflation_rate) ** year_number
        total += (price * curr_inflation) * 12
    return total


def roi(zestimate, inflation_rate, mortgage_rate, num_years, rental_price, down_payment_percent):
    down_payment = zestimate * down_payment_percent
    P = zestimate * (1 - down_payment_percent)

    incurred_cost = calculate_monthly(
        P, mortgage_rate, num_years) * 12 * num_years + down_payment

    income = airbnb_income(price=rental_price,
                           inflation_rate=inflation_rate,
                           num_years=num_years)

    return (income - incurred_cost) / incurred_cost


# -

# ### Model 1
#
# - Simpler causal model without exogenous variables

# +
Neigh_alias = ['Neigh_North', 'Neigh_South']
BPB_alias = [2, 3, 4, 5, 6]

north_prob = float(len(df[df['Neigh'] == 'Neigh_North']))/len(df)
bpb_prob = [(df.BPB == a).count()/len(df) for a in BPB_alias]


Neigh_prob = torch.tensor([north_prob, 1 - north_prob])
BPB_prob = torch.tensor(bpb_prob)


def model():
    Neigh = pyro.sample("Neigh", dist.Categorical(probs=Neigh_prob))
    BPB = pyro.sample("BPB", dist.Categorical(probs=BPB_prob))

    bpb_converted = BPB_alias[BPB]
    neigh_converted = Neigh_alias[Neigh]
    rent_pred = rent_model.predict(pd.DataFrame({'BPB': [bpb_converted]}))[0]
    zest_pred = zest_model.predict(pd.DataFrame(
        {'BPB': [bpb_converted], 'Neigh': [neigh_converted]}))[0]

    Rent = pyro.sample("Rent", dist.Delta(torch.tensor(rent_pred)))
    Zest = pyro.sample("Zest", dist.Delta(torch.tensor(zest_pred)))

    ROI = pyro.sample("ROI", dist.Delta(roi(zestimate=Zest,
                                            rental_price=Rent,
                                            inflation_rate=INFLATION_RATE,
                                            mortgage_rate=MORTGAGE_RATE,
                                            num_years=NUM_YEARS,
                                            down_payment_percent=DOWN_PAYMENT)))

    return {'Neigh': Neigh, 'BPB': BPB, 'Rent': Rent, 'Zest': Zest, 'ROI': ROI}


# -

# ### Model 2 - SCM

# \begin{align}
# Nb &= p(BPB=2)=p_2, p(BPB=3)=p_3, p(BPB=4)=p_4,p(BPB=5)=p_5, p(BPB=6)=p_6 \\
# Nr &= \mathcal{N}(0, \text{rent model standard error}) \\
# Nz &= \mathcal{N}(0, \text{zest model standard error}) \\
# Nn &= Bernoulli(p_{south})\\\\
# Neigh &= Nn \\
# BPB &= Nb \\
# Rent &= bpb_3 \cdot a + bpb_4 \cdot b + bpb_5 \cdot d + bpb_6 \cdot d + e +Nrent\\
# Zest &= bpb_3 \cdot a + bpb_4 \cdot b + bpb_5 \cdot d + bpb_6 \cdot d + NeighSouth \cdot e + f + Nz \\
# ROI &= \text{ROIfunc}(Rent, Zest)
# \end{align}

# +
Neigh_alias = ['Neigh_North', 'Neigh_South']
BPB_alias = [2, 3, 4, 5, 6]

south_prob = float(len(df[df['Neigh'] == 'Neigh_South']))/len(df)
bpb_prob = [(df.BPB == a).count()/len(df) for a in BPB_alias]


# Neigh_prob = torch.tensor([north_prob, 1- north_prob])
BPB_prob = torch.tensor(bpb_prob)

exogenous_dists = {
    'Nn': dist.Bernoulli(torch.tensor(south_prob)),
    'Nb': dist.Categorical(BPB_prob),
    'Nr': dist.Normal(torch.tensor(0.), torch.tensor(rent_model_se)),
    'Nz': dist.Normal(torch.tensor(0.), torch.tensor(zest_model_se))
}


def scm(exogenous_dists):
    Nn = pyro.sample("Nn", exogenous_dists['Nn'])
    Nb = pyro.sample("Nb", exogenous_dists['Nb'])
    Nr = pyro.sample("Nr", exogenous_dists['Nr'])
    Nz = pyro.sample("Nz", exogenous_dists['Nz'])

    z_bounds = torch.tensor(zest_model_se)
    r_bounds = torch.tensor(rent_model_se)

    # Stay within 1 SD of 0
    Nz = min(z_bounds, Nz)
    Nz = max(-2 * z_bounds, Nz)

    Nr = min(r_bounds, Nr)
    Nr = max(-2 * r_bounds, Nr)

    Neigh = pyro.sample("Neigh", dist.Delta(Nn))
    BPB = pyro.sample("BPB", dist.Delta(Nb))

    bpb_converted = BPB_alias[BPB.int()]
    neigh_converted = Neigh_alias[Neigh.int()]

    rent_pred = rent_model.predict(pd.DataFrame({'BPB': [bpb_converted]}))[0]
    zest_pred = zest_model.predict(pd.DataFrame(
        {'BPB': [bpb_converted], 'Neigh': [neigh_converted]}))[0]

    Rent = pyro.sample("Rent", dist.Delta(torch.tensor(rent_pred) + Nr))
    Zest = pyro.sample("Zest", dist.Delta(torch.tensor(zest_pred) + Nz))

    ROI = pyro.sample("ROI", dist.Delta(roi(zestimate=Zest,
                                            rental_price=Rent,
                                            inflation_rate=INFLATION_RATE,
                                            mortgage_rate=MORTGAGE_RATE,
                                            num_years=NUM_YEARS,
                                            down_payment_percent=DOWN_PAYMENT)))

    return {'Neigh': Neigh, 'BPB': BPB, 'Rent': Rent, 'Zest': Zest, 'ROI': ROI}


# -

# # Model Analysis

# ### Causal Effect on ROI of Changing Neighborhoods - Should we choose North or South?
#
# $$
# P(ROI | do(Neigh=South)) - P(ROI|do(Neigh=North))
# $$

# +
south_model = pyro.do(
    model, data={'Neigh': torch.tensor(Neigh_alias.index('Neigh_South'))})
north_model = pyro.do(
    model, data={'Neigh': torch.tensor(Neigh_alias.index('Neigh_North'))})

south_scm = pyro.do(
    scm, data={'Neigh': torch.tensor(Neigh_alias.index('Neigh_South'))})
north_scm = pyro.do(
    scm, data={'Neigh': torch.tensor(Neigh_alias.index('Neigh_North'))})
# -

# Causal Effect for Model 2(SCM)

# +
scm_south_roi_samples = [south_scm(exogenous_dists)[
    'ROI'].item() for _ in range(2000)]
scm_north_roi_samples = [north_scm(exogenous_dists)[
    'ROI'].item() for _ in range(2000)]

scm_causal_effect_neigh = mean(
    scm_south_roi_samples) - mean(scm_north_roi_samples)
scm_causal_effect_neigh
# -

# Showing the problematic large variance in the samples

mean(scm_south_roi_samples), np.std(scm_south_roi_samples), min(
    scm_south_roi_samples), max(scm_south_roi_samples)

# Causal Effect for Model 1

# +
south_roi_samples = [south_model()['ROI'].item() for _ in range(2000)]
north_roi_samples = [north_model()['ROI'].item() for _ in range(2000)]

causal_effect_neigh = mean(south_roi_samples) - mean(north_roi_samples)
causal_effect_neigh
# -

# ### Results:
# - SCM Model has too much variance and so it makes it difficult to get good samples
# - Model 1 shows a 1.3% Increase In ROI from Investing in property in the South as opposed to the North
#     - Why?
#         - South includes: San Jose, Cupertino and other properties a bit further from SF
#         - North includes: Mountain View, Palo Alto - properties in the heart of Silicon Valley and closer to SF

# ### Causal Effect on ROI of Different BPBs

bpb_roi_samples = []
for i in range(len(BPB_alias)):
    bpb = BPB_alias[i]
    bpb_model = pyro.do(model, data={'BPB': torch.tensor(i)})
    bpb_roi_samples.append([bpb_model()['ROI'].item() for _ in range(2000)])

for i in range(len(BPB_alias)):
    bpb = BPB_alias[i]
    print(bpb, mean(bpb_roi_samples[i]))

# ### What about BPB=5 vs BPB=2 ?
# $P(ROI | do(BPB=5)) - P(ROI | do(BPB=2))$

mean(bpb_roi_samples[3]) - mean(bpb_roi_samples[0])

# ### Results 
# - 3 Beds 2 baths is the best configuration for maximizing ROI
# - You can increase ROI by 1.5% by doing 3 Beds 2 Baths over a 1 Bed 1 Bath

# ### Given we want a certain Bed and Baths, where should we buy?
#
# Example for 2 bed 1 bath:
# $$
# P(ROI | do(Neigh= North), BPB=3) - P(ROI | do(Neigh=South), BPB=3)
# $$

for bpb_index in range(len(BPB_alias)):
    bpb = BPB_alias[bpb_index]
    roi_means = []
    cond_model = pyro.condition(model, {"BPB": bpb_index})

    for neigh_index in range(len(Neigh_prob)):
        curr_model = pyro.do(cond_model, {'Neigh': neigh_index})
        roi_mean = curr_model()['ROI'].item()
        roi_means.append(roi_mean)

    print("𝑃(𝑅𝑂𝐼|𝑑𝑜(𝑁𝑒𝑖𝑔ℎ=𝑁𝑜𝑟𝑡ℎ),𝐵𝑃𝐵=%d)−𝑃(𝑅𝑂𝐼|𝑑𝑜(𝑁𝑒𝑖𝑔ℎ=𝑆𝑜𝑢𝑡ℎ),𝐵𝑃𝐵=%d) = %.3f - %.3f = %.3f" % (bpb, bpb,
                                                                                                roi_means[0],
                                                                                                roi_means[1],
                                                                                                roi_means[0] -
                                                                                                roi_means[1]
                                                                                                ))

# ### Results
# - No matter which BPB you are looking at, it is always better to pick a property in the South rather than the North

# ### Given we want to buy in a certain area, what beds and baths should we buy?
#
# Example for the North Neighborhood
# $$
# P(ROI | do(BPB= x), Neigh= North) - P(ROI | do(BPB=y), Neigh= North)
# $$

for neigh_index in range(len(Neigh_prob)):
    roi_means = []
    cond_model = pyro.condition(model, {"Neigh": neigh_index})
    neigh = Neigh_alias[neigh_index]
    print(neigh)

    for bpb_index in range(len(BPB_prob)):
        bpb = BPB_alias[bpb_index]
        curr_model = pyro.do(cond_model, {'BPB': bpb_index})
        roi_mean = curr_model()['ROI'].item()
        print("𝑃(𝑅𝑂𝐼|𝑑𝑜(𝐵𝑃𝐵=%s),𝑁𝑒𝑖𝑔ℎ=%s) = %.3f" % (bpb, neigh, roi_mean))

# ### Results
# - If you are interested in purchasing a property in the North, it should be a 3 bed 2 bath. If you are interested in purchasing a property in the South it should also be a 3 bed 2 bath.

# # Property Examples
# - here are the properties with the highest ROIs
# - We can see that these properties follow the effects we have seen in our model
#     - 5 BPB, all are South Neighborhood

df.sort_values(["ROI"], ascending=False)[:5]

# - [Property 1](https://www.zillow.com/homedetails/435-Milan-Dr-San-Jose-CA-95134/2143393753_zpid/)
# - [Property 2](https://www.zillow.com/homedetails/3215-Lamond-Ct-San-Jose-CA-95148/19793643_zpid/)
# - [Property 5](https://www.zillow.com/homedetails/2409-Venturi-Pl-UNIT-2-San-Jose-CA-95132/124743156_zpid/)

# # Next Steps
#
# 1. **Get more data**
# 2. Consider adding new nodes to the DAG
#     1. Individual Neighborhoods as opposed to North and South
#     2. Host and ratings
#     3. Distance to transit
# 3. Reconsider new edges with our new data
# 3. Test more complex functions to fit Rent and Zestimate
#     1. Currently using factor effects models
# 4. Attempt to use SCM again
#     1. Our models have high Std Error so using SCM is difficult
#     2. Increasing the amount of data would also decrease Std. Error and help give better results from SCM

StarcoderdataJupyterScriptsDedupFiltered