Buenos_Aires_Price_Prediction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9c4f669e-a914-47b9-9767-c1e9419125ba",
   "metadata": {},
   "source": [
    "<h1 style=\"text-align: center; font-family: 'Times New Roman', serif; font-weight: bold; font-size: 60px;\">Predicting Apartments Prices in Buenos Aires, Argentina</h1>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6773f96",
   "metadata": {},
   "source": [
    "<div style=\"width: 90%; font-family: Times New Roman, serif; text-justify: inter-word; margin: 0 auto; font-size: 20px;\">\n",
    "\n",
    "## **Introduction**\n",
    "### **Motivation**\n",
    "The real estate market in ***Buenos Aires*** is influenced by various factors such as economic conditions, population growth, and urban development. Predicting house prices accurately is essential for various stakeholders, including buyers, sellers, real estate agents, and policy makers. Accurate price predictions can help in making informed decisions, optimizing investments, and planning sustainable urban development.\n",
    "\n",
    "\n",
    "### **Objective**\n",
    "We aim to develop a reliable model to predict apartment prices in Buenos Aires using available data, despite the absence of temporal indicators. \n",
    "\n",
    "## **Table of Contents**\n",
    "1. [Introduction](#introduction)\n",
    "2. [Problem Statement](#problem-statement)\n",
    "3. [Data Description](#data-description)\n",
    "4. [Methodology](#methodology)\n",
    "5. [Exploratory Data Analysis (EDA)](#exploratory-data-analysis-eda)\n",
    "6. [Data preparation (preprocessing)](#data-preparation-preprocessing)\n",
    "7. [Modeling](#modeling)\n",
    "8. [Reverse Heteroskedasticity](#reverse-heteroskedasticity)\n",
    "9. [Final Model: Gradient Boosting Regressor](#final-model-gradient-boosting-regressor)\n",
    "10. [Store Model](#store-model)\n",
    "11. [Deployment](#deployment)\t\n",
    "12. [Conclusion](#conclusion)\n",
    "13. [References](#references)\n",
    "14. [About the Author](#about-the-author)\n",
    "\n",
    "## **Problem Statement**\n",
    "\n",
    "Identify the significant features that accurately predict apartment prices in Buenos Aires for units costing less than $400,000 USD using a suitable regression model. Aim to achieve a mean absolute error (MAE) of less than 50% compared to the baseline model, using a dataset of 4,000+ apartment listings.\n",
    "\n",
    "## **Data Description**\n",
    "\n",
    "### **Data Sources**\n",
    "The dataset for this project is sourced from websites that provide real estate services in Buenos Aires. It includes 12,000+ house listings from Buenos Aires, containing various features relevant to housing prices. The data was scraped to ensure comprehensive and up-to-date information, capturing the unique characteristics of the local real estate market.\n",
    "\n",
    "### **Data Overview**\n",
    "\n",
    "The dataset consists of over 43,000 listings with 14 features. Below is a description of each feature:\n",
    "\n",
    "- **operation**: Type of operation (e.g., sale, rent).\n",
    "- **property_type**: Type of property (e.g., house, apartment, PH (Propiedad Horizontal), Store).\n",
    "- **place_with_parent_names**: The location of the property, including parent regions for context.\n",
    "- **lat-lon**: Latitude and longitude coordinates of the property.\n",
    "- **price**: The listed price of the property.\n",
    "- **currency**: The currency of the listed price (e.g., USD, ARS).\n",
    "- **price_aprox_local_currency**: Approximate price in the local currency.\n",
    "- **price_aprox_usd**: Approximate price in USD.\n",
    "- **surface_total_in_m2**: Total surface area of the property in square meters.\n",
    "- **surface_covered_in_m2**: Covered surface area of the property in square meters.\n",
    "- **price_usd_per_m2**: Price per square meter in USD.\n",
    "- **price_per_m2**: Price per square meter in the listed currency.\n",
    "- **floor**: The floor number of the property (if applicable).\n",
    "- **rooms**: The number of rooms in the property.\n",
    "- **expenses**: Expenses associated with the property.\n",
    "- **properati_url**: URL of the property listing on the Properati website.\n",
    "\n",
    "\n",
    "## **Methodology**\n",
    "\n",
    "This project follows a prescriptive methodology to predict apartment prices in Buenos Aires, starting with defining the objective of identifying significant features using a suitable regression model. Key metrics were set, including achieving a mean absolute error (MAE) of less than 50% compared to a baseline model. \n",
    "\n",
    "The data, scraped from real estate websites, was cleaned by handling missing values, converting data types, removing duplicates, and normalizing features. Exploratory Data Analysis (EDA) was conducted to understand feature distributions and relationships.\n",
    "\n",
    "Feature engineering was used at each step to identify the most significant features in our model. \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae9d0cdb",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; \">\n",
    "\n",
    "## **Exploratory Data Analysis (EDA)**\n",
    "\n",
    "We thoroughly examine the data to determine the preprocessing steps that come after.\n",
    "\n",
    "### **Import Libraries**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0d1c9e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import plotly.express as px\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import missingno as msno\n",
    "import scipy.stats as sps\n",
    "\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.model_selection import train_test_split\n",
    "import statsmodels.api as sm\n",
    "from sklearn.metrics import mean_absolute_error\n",
    "\n",
    "warnings.simplefilter(action=\"ignore\", category=FutureWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2500129",
   "metadata": {},
   "source": [
    "\n",
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "### **Import Data**\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0cd36385",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read and Convert the CSV file into a Pandas DataFrame\n",
    "df_buenos_aires = pd.read_csv('data/buenos-aires-real-estate.csv', low_memory=False)\t"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "daf36813",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "### **Investigate Data**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49a28842",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_buenos_aires.info()\t"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf121fab",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "The output of the method *info()* provides many information. For instance, we clearly observe that the data is not complete and missing values are prevalent is some columns. In addition, both columns `expenses` and `lat-lon` should not be of type *object*. \n",
    "\n",
    "#### **Missing Values**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2464e369",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Percentage of missing values by column\n",
    "nan_columns_percentage = (df_buenos_aires.isna().sum(axis=0) * 100 / len(df_buenos_aires)).sort_values(ascending=True)\n",
    "\n",
    "plt.figure(figsize=(10, 6))\n",
    "sns.barplot(x=nan_columns_percentage, y=nan_columns_percentage.index, orient='h', \n",
    "            color='royalblue', alpha=0.8, edgecolor='black')\n",
    "plt.title('Percentage of missing values by column')\n",
    "plt.xlabel('Percentage (%)')\n",
    "plt.ylabel('Column')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5b2a288",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "Features like `expenses` and `floor` have more than 50% of missing values, which informs our decision to just drop them, because no matter what filling strategy we use, we are prone to influence our predictive model and therefore affect accuracy. \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af924de4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot missing values with color royalblue\n",
    "msno.matrix(df_buenos_aires, color=(0.3, 0.4, 0.8));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7494ebe8",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "The missing value matrix provide us with many information in a clear and concise way. We can observe that whenever `price` is missing, `currency`, `price_aprox_local_currency`, and `price_aprox_usd` are also missing. Which suggests an association between them.\n",
    "\n",
    "The `expenses` column is the one with the most missing values, followed by the `floor` and `rooms` and then `price_usd_per_m2`.\n",
    "\n",
    "#### **Describe Data**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c4c302d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Describe the dataset\n",
    "df_buenos_aires.describe(include='all').T"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a73c2649",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "The method `describe()` provide a descriptive analysis of our data. We can observe that `price`, `price_aprox_local_currency`, `price_aprox_usd`, `surface_total_in_m2`, `surface_covered_in_m2`, and `price_usd_per_m2` show a minimum value of **0.0** which is not likely in our case. On the other hand, having a **987897.0 $m^2$** in the `surface_covered_in_m2` is exaggerated and could be a sign for an outlier.\n",
    "\n",
    "These columns may have outliers we should deal with later.\n",
    "\n",
    "#### **Histograms and boxplots**\n",
    "In this analysis, we focus on the surface covered feature, which represents the total indoor living area of a property, including all enclosed spaces such as rooms, hallways, bathrooms, and kitchens. \n",
    "\n",
    "This metric is particularly important from the perspective of homebuyer considerations, as it directly impacts the usability and comfort of the living space. Homebuyers often prioritize the amount of indoor space available for daily activities, family gatherings, and personal comfort. \n",
    "\n",
    "By analyzing the surface covered, we can provide more accurate insights into the value and livability of apartments, helping potential buyers make informed decisions based on the actual indoor area they will be utilizing. This focus ensures that our model aligns with the key factors that homebuyers consider when evaluating properties.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a538eb38",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot distribution of apartment sizes \n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))\n",
    "\n",
    "# Distribution of surface covered\n",
    "df_buenos_aires.surface_covered_in_m2.plot(kind='hist', color='royalblue', alpha=0.8, edgecolor='black', ax=ax1)\n",
    "\n",
    "# Boxplot of surface covered \n",
    "#df_buenos_aires.surface_covered_in_m2.plot(kind='box', color='royalblue', vert=False, rot=90, ax=ax2)\t\n",
    "sns.boxenplot(data=df_buenos_aires, x='surface_covered_in_m2', color='royalblue', ax=ax2, scale='linear')\n",
    "\n",
    "ax1.set_xlabel('Area [sq meters]')\n",
    "ax2.set_xlabel('Area [sq meters]')\n",
    "fig.suptitle('Distribution of Apartment Sizes (Surface Covered)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0a2201a",
   "metadata": {},
   "source": [
    "<div style=\" font-family: 'Times New Roman', serif;\">\n",
    "Apparently, most of our listings have an Area of no more than 100,000 $m^2$. The boxplot confirms the existence of outliers.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75c03ca8",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_buenos_aires.surface_covered_in_m2.describe().to_frame()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56960200",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "Based on the Histogram the boxplot and the descibe method for the `surface_covered_in_m2` column, we conclude that there are outliers in the data. The strategy chosen in this analysis is to drop them and we will use the *quantile strategy*, the quantiles are to be set after exploring the data.\n",
    "\n",
    "### **Quantile Strategy**\n",
    "\n",
    "#### **Quantile Categorization**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fcec4605",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_quantile_category(df, column, quantile_range):\n",
    "    \"\"\"\n",
    "    Categorize the data into quantile categories and return a DataFrame\n",
    "    with the quantile category, count, min, and max per quantile category.\n",
    "    \n",
    "    Args: \n",
    "        df: Pandas DataFrame\n",
    "        column: Column to categorize\n",
    "        quantiles: List of quantiles\n",
    "    \n",
    "    Returns:\n",
    "        agg_df: Pandas DataFrame with quantile categories, count, min, and max\n",
    "    \"\"\"\n",
    "\n",
    "    quantiles = df[column].quantile(quantile_range)\n",
    "\n",
    "    # Categorize the data\n",
    "    n = len(quantiles.unique())\n",
    "    labels = quantile_range[1:n]\n",
    "\n",
    "    categories = pd.cut(df[column], bins=quantiles.values, labels=labels, include_lowest=True, retbins=False, right=True, duplicates='drop')\n",
    "    categories = categories.rename('quantile_category')\n",
    "    df_quantile = pd.concat([df, categories], axis=1)\n",
    "    \n",
    "    # Calculate the count, min, and max for each quantile category\n",
    "    agg_df = df_quantile.groupby('quantile_category')[column].agg(['count', 'min', 'max']).reset_index()\n",
    "    \n",
    "    return agg_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ec82889",
   "metadata": {},
   "outputs": [],
   "source": [
    "quantile_range = np.linspace(0, 1, 30).round(2) # 29 quantiles\n",
    "df_quantiles_min_max = get_quantile_category(df_buenos_aires, 'surface_covered_in_m2', quantile_range)\n",
    "df_quantiles_min_max"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b31fc29c",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "Based on the results of the quantile table, we can see that the **10% of the smallest values** are **<= 34** sq meters and the *2*% of the largest values** are **>= 500** sq meters.\n",
    "\n",
    "We can observe clearly that the sudden jump after the **98th percentile** is big compared to the rest of the data. That's why we decided to drop these outliers.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e589783",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5), sharey=True)\n",
    "\n",
    "# Plot distribution of apartment sizes greater than 450 sq meters \n",
    "(df_buenos_aires[df_buenos_aires.surface_covered_in_m2 >= 450]\n",
    " .surface_covered_in_m2\n",
    " .plot(kind='hist', edgecolor='black', color='royalblue', bins=500, alpha=0.8, ax=ax1, xticks=range(0, 12_900, 600), xlim=(0, 12_900), rot=45))\n",
    "\n",
    "# Plot distribution of apartment sizes less than 20 sq meters\n",
    "(df_buenos_aires[df_buenos_aires.surface_covered_in_m2 <= 30]\n",
    " .surface_covered_in_m2\n",
    " .plot(kind='hist', edgecolor='black', color='royalblue', alpha=0.8, ax=ax2, xticks=range(0, 31, 2)))\n",
    "\n",
    "ax1.set_xlabel('Area [sq meters]')\n",
    "ax2.set_xlabel('Area [sq meters]')\n",
    "\n",
    "plt.suptitle('Distribution of Apartment Sizes (Surface Covered)', y=1.05, size=25)\n",
    "\n",
    "ax1.set_title('\\nGreater than 450 sq meters', size=18)\n",
    "ax2.set_title('\\nLess than 20 sq meters', size=18)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd170ff0",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "The histogram on the right tells us that more than 900 observations have an area less then 3 sq meters, which is not expected. Additionally, we can see that there are about 400 observations with an area between 18 and 21 sq meters. \n",
    "\n",
    "Similarly, the histogram on the left tells us that there are more than 1,000 observations with an area greater than 450 and less than 2,400 sq meters. \n",
    "The histograms confirm what we expected, therefore we will proceed with the chosen quantile strategy.\n",
    "\n",
    "**NOTE:** The scale of the x-axis in the left histogram starts at 450 sq meters and ends at 125,000 sq meters.\n",
    "\n",
    "#### **Remove Outliers (quantile strategy)**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61a00fe6",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def remove_outliers(df, column, q_a, q_b):\n",
    "    \"\"\"\n",
    "    Remove outliers using the percentiles\n",
    "    Args: \n",
    "        df: Pandas DataFrame\n",
    "        column: Column to remove outliers   \n",
    "    Returns:\n",
    "        df: Pandas DataFrame with outliers removed\n",
    "    \"\"\"\n",
    "    qa_ = df[column].quantile(q_a)\n",
    "    qb_ = df[column].quantile(q_b)\n",
    "    return df[(df[column].between(qa_, qb_, inclusive='both'))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89751462",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate quantiles;\n",
    "q1, q2 = \tdf_buenos_aires.surface_covered_in_m2.quantile([0.10, 0.98])\n",
    "print(f'The 10th and 98th percentile are {q1: .2f} and {q2: .2f} respectively')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c69057e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot distribution of apartment sizes after removing outliers\n",
    "remove_outliers(df_buenos_aires, 'surface_covered_in_m2', 0.1, 0.9).surface_covered_in_m2.describe().to_frame() # type: ignore"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "511ff1b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot distribution of apartment sizes after removing outliers\t\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))\n",
    "\n",
    "# Distribution of surface covered\n",
    "remove_outliers(df_buenos_aires, 'surface_covered_in_m2', 0.1, 0.9).surface_covered_in_m2.plot(kind='hist', color='royalblue', alpha=0.8, edgecolor='black', ax=ax1)\n",
    "\n",
    "# Boxen plot of surface covered, you can choose the width method to be either `linear` or `exponential`\n",
    "sns.boxenplot(data=remove_outliers(df_buenos_aires, 'surface_covered_in_m2', 0.1, 0.9), x='surface_covered_in_m2', color='royalblue', ax=ax2, width_method='exponential')\n",
    "\n",
    "ax1.set_xlabel('Area [sq meters]')\n",
    "ax2.set_xlabel('Area [sq meters]')\n",
    "plt.suptitle('Distribution of Apartment Sizes (Surface Covered) \\nAfter Removing Outliers')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffb3e79f",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "Now the distribution of apartment sizes looks much more natural with the outliers removed.\n",
    "\n",
    "The histogram is right skewed, which is expected. As less and less apartments have larger surface area. \n",
    "\n",
    "The `boxenplot` in the right, shows the distribution of quantiles for the `surface_covered_in_m2` column. Each successive level outward contains half of the remaining data. So the first two sections out from the centerline contain 50% of the data. After that, the next two sections contain 25% of the data. This continues until we are at the outlier level. Each level out is shaded lighter.\n",
    "\n",
    "We observe that the right side of each level contains more data than the left side. This further confirms what we expected.\n",
    "\n",
    "Following the same approach, we can remove outliers from the `surface_total_in_m2` column. We will use the same quantile strategy, since the two features are closely related.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "527e8b5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "get_quantile_category(df_buenos_aires, 'surface_total_in_m2', quantile_range)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "278bd2b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "remove_outliers(df_buenos_aires, 'surface_total_in_m2', 0.2, 0.9).surface_total_in_m2.describe().to_frame() # type: ignore"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4df98fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot distribution of apartment sizes (Surface Total) after removing outliers\t\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))\n",
    "\n",
    "# Distribution of surface covered\n",
    "remove_outliers(df_buenos_aires, 'surface_total_in_m2', 0.2, 0.9).surface_total_in_m2.plot(kind='hist', color='royalblue', alpha=0.8, edgecolor='black', ax=ax1)\n",
    "\n",
    "# Boxen plot of surface covered \n",
    "sns.boxenplot(data=remove_outliers(df_buenos_aires, 'surface_total_in_m2', 0.2, 0.9), x='surface_total_in_m2', color='royalblue', ax=ax2, width_method='exponential')\n",
    "\n",
    "ax1.set_xlabel('Area [sq meters]')\n",
    "ax2.set_xlabel('Area [sq meters]')\n",
    "plt.suptitle('Distribution of Apartment Sizes (Surface Total) \\nAfter Removing Outliers')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c28056c",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "Again, following the same approach, we can remove outliers from the `rooms` column. We will use the same quantile strategy. \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "225477ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "get_quantile_category(df_buenos_aires, 'rooms', quantile_range)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2644c2a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_buenos_aires.rooms.describe().to_frame()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf39b309",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_buenos_aires.rooms.plot(kind='hist', edgecolor='black', color='royalblue', bins=30, alpha=0.8, xticks=range(0, 32, 2), xlim=(0, 33))\t\n",
    "plt.xlabel('Number of Rooms')\n",
    "plt.title('Distribution of Number of Rooms')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "850e97d7",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "It looks like most of our apartments have between 1 and 7 rooms. Only few apartments have more than 7 rooms, and again this was expected.\n",
    "\n",
    "Moving on, we will explore the distribution of `price_aprox_usd` using the approach as before.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c6c890f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot distribution of apartment prices\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))\n",
    "\n",
    "# Distribution of surface covered\n",
    "df_buenos_aires.price_aprox_usd.plot(kind='hist', color='royalblue', alpha=0.8, edgecolor='black', ax=ax1)\n",
    "\n",
    "# Boxen plot of surface covered \n",
    "sns.boxenplot(data=df_buenos_aires, x='price_aprox_usd', color='royalblue', ax=ax2, scale='exponential')\n",
    "\n",
    "ax1.set_xlabel('Price [USD]')\n",
    "ax2.set_xlabel('Price [USD]')\n",
    "plt.suptitle('Distribution of Apartment Prices')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "653cd1f9",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "We can see that most of our listings have an approximate price of less than $800,000. The `boxenplot` confirms the existence of outliers.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d128d257",
   "metadata": {},
   "outputs": [],
   "source": [
    "price_quantile_range = np.linspace(0, 1, 90).round(2)\n",
    "get_quantile_category(df_buenos_aires, 'price_aprox_usd', price_quantile_range)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "152781a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "remove_outliers(df_buenos_aires, 'price_aprox_usd', 0.1, 0.9).price_aprox_usd.describe().to_frame() # type: ignore"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfa777b3",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "We will use the 1st and the 97th quantile to remove outliers. There's a sudden jump in the prices of appartments after the 97th quantile, if left out it would affect the distribution of prices and thus our predictive model.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f2084ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot distribution of apartment prices after removing outliers\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))\n",
    "\n",
    "# Distribution of surface covered\n",
    "remove_outliers(df_buenos_aires, 'price_aprox_usd', 0.1, 0.9).price_aprox_usd.plot(kind='hist', color='royalblue', alpha=0.8, edgecolor='black', ax=ax1)\n",
    "\n",
    "# Boxen plot of surface covered \n",
    "sns.boxenplot(data=remove_outliers(df_buenos_aires, 'price_aprox_usd', 0.1, 0.9), x='price_aprox_usd', color='royalblue', ax=ax2, scale='exponential')\n",
    "\n",
    "ax1.set_xlabel('Area [sq meters]')\n",
    "ax2.set_xlabel('Area [sq meters]')\n",
    "plt.suptitle('Distribution of Apartment Prices \\nAfter Removing Outliers')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d8b3b60",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "### **Data Exploration**\n",
    "#### **Analysis of correlation**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a546fbac",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select only numerical columns\n",
    "numerical_columns = df_buenos_aires.select_dtypes(include='number').columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1adcb53f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Correlation matrix between numerical columns\n",
    "corr_df = df_buenos_aires[numerical_columns].corr()  # type: ignore\t\n",
    "\n",
    "# Mask the upper triangle of the correlation matrix\n",
    "mask_corr = np.triu(np.ones_like(corr_df, dtype=bool))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0062d8b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize = (12, 8))\n",
    "sns.heatmap(corr_df, annot = True, cmap = \"Spectral\", linewidths = 2, linecolor = \"white\", fmt='.2f', mask=mask_corr)\n",
    "plt.title(\"Correlation Matrix\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "684998ea",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "The first thing we notice from this correlation matrix is that `price_aprox_usd` and `price_aprox_local_currency` are perfectly correlated, which is expected because we can get `price_aprox_usd` from `price_aprox_local_currency` and exchange rate. \n",
    "\n",
    "The second thing we notice is that `price_aprox_local_currency` and `price_aprox_usd` are somewhat moderately correlated with `rooms`. This suggests that the higher the number of rooms, the higher the price.\n",
    "\n",
    "The third thing we notice is that `price_aprox_local_currency` and `price_aprox_usd` are moderately correlated with `price`. \n",
    "\n",
    "The fourth thing we notice is that `price_usd_per_m2` and `rooms` have a weakly negative correlation. Which is not expected. \n",
    "`price`, `price_aprox_local_currency`, and `price_aprox_usd` can present high **co-linearity**, and thus affect the accuracy of our predictive model. We'll have to only keep one of them, we choose `price_aprox_usd` for our model.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70961806",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot distribution of property types\n",
    "plt.figure(figsize=(8, 8))\n",
    "df_buenos_aires.property_type.value_counts().plot(kind='pie', autopct='%1.1f%%', shadow=True, explode=[0, 0, 0.2, 0.3])\n",
    "plt.ylabel('')\n",
    "plt.title('Distribution of Property Types')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b9c1143",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "We can see that most of our listings are of `apartment` type and `house` type representing around 90% of our data.  \n",
    "\n",
    "We should consider filtering our data to only include one type of property in our model to avoid inaccuracy.\n",
    "\n",
    "#### **Investigate Data time period**\n",
    "\n",
    "We will try to investigate the data period of the dataset. We will use the `price_aprox_local_currency` and `price_aprox_usd` to find the exchange rate and then **reverse-engineer** that information with the help of public exchange rates to find precisely when the data was collected.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90d28114",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate exchange rate between local currency and USD\n",
    "(df_buenos_aires.price_aprox_local_currency / df_buenos_aires.price_aprox_usd).round(5).unique()[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "184fd738",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "We found that the exchange rate was 15.1624, and by referencing the public exchange rates, we conclude that the data was collected between **Sepetember** and\t**October 2016**. \n",
    "\n",
    "#### **Investigate currency**\n",
    "\n",
    "We will drop all the NaN values based on the `price` column.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "154a6d4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop NaN Values.\n",
    "df_buenos_aires = df_buenos_aires[~df_buenos_aires.price.isna()]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93a16a6f",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "We will drop observations that have `price` equal to `0.0`.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68a0b811",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop rows where price is equal to 0\n",
    "df_buenos_aires = df_buenos_aires[df_buenos_aires.price != 0.0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "397a44d9",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "Now will invistigate prices that have a `USD` currency and have a different value from the one shown in the column `price_aprox_usd`\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74155e8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check if price_aprox_usd is not equal to price when currency is USD\n",
    "df_buenos_aires[(df_buenos_aires.price != df_buenos_aires.price_aprox_usd) & (df_buenos_aires.currency == 'USD')].shape[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "036108b3",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "#### **Investigate price_usd_per_m2 and price_aprox_usd**\n",
    "\n",
    "We will investigate the relationship between `price_usd_per_m2`, and `price_aprox_usd`.\n",
    "\n",
    "It seems that the `price_usd_per_m2` can be derived from the `price_aprox_usd` and `surface_covered_in_m2`.\n",
    "Similarly, the `price_per_m2` can be derived from the `price_aprox_local_currency` and `surface_total_in_m2`.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c034b8da",
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.options.display.float_format = '{:.2f}'.format\n",
    "\n",
    "# we can get price_usd_per_m2 by dividing price_aprox_usd by surface_total_in_m2\n",
    "(df_buenos_aires.price_usd_per_m2\n",
    "- df_buenos_aires.price_aprox_usd\n",
    "/ df_buenos_aires.surface_total_in_m2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ea41be5",
   "metadata": {},
   "outputs": [],
   "source": [
    "(df_buenos_aires.price_usd_per_m2\n",
    "- df_buenos_aires.price_aprox_usd\n",
    "/ df_buenos_aires.surface_total_in_m2).describe().to_frame()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a912075",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "#### **Investigate place_with_parent_names**\n",
    "\n",
    "In this part, we will extract the country, province, municipality, and neighborhood from the `place_with_parent_names` column. \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0051742b",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_buenos_aires.place_with_parent_names.str.split('|', expand=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4feb5bf5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split the `place_with_parent_names` column into multiple columns\n",
    "df_localities = df_buenos_aires.place_with_parent_names.str.split('|', expand=True)\n",
    "\n",
    "# Drop columns 0, 5, 6\n",
    "df_localities.drop(columns=[0, 5, 6], inplace=True)\n",
    "df_localities.rename(columns={1:'country', 2:'province', 3:'municipality', 4:'neighborhood'}, inplace=True)\n",
    "\n",
    "# Concatenate `df_buenos_aires` and `df_localities`\n",
    "df_buenos_aires = pd.concat([df_buenos_aires, df_localities], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7e8d242",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Count listings per province\t\n",
    "listings_per_province = df_buenos_aires.groupby(by='province').province.value_counts().sort_values(ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c935c0fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 8))\n",
    "\n",
    "((listings_per_province * 100 / len(df_buenos_aires))\n",
    ".plot.pie(explode=[0.01, 0.01, 0.05, 0.1], autopct='%1.1f%%', shadow=True,\n",
    "            startangle=45, labeldistance=None,\n",
    "            colors=['royalblue', 'crimson', 'limegreen', 'orange'])\n",
    "            )   # type: ignore\n",
    "\n",
    "plt.ylabel('')\n",
    "plt.legend(listings_per_province.index, loc=\"lower left\")\n",
    "plt.title('Percentage of listings per province')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c501df5",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "We can see that province `Bs.As. G.B.A. Zona Norte` has the highest listings with 38%. Followed by `Capital Federal` with 30% and `Bs.As. G.B.A. Zona Oeste` with 21%. Finally, `Bs.As. G.B.A. Zona Sur` has the lowest listings with 9%.\n",
    "A quick research tells us that the province with the highest listings (`Bs.As. G.B.A. Zona Norte`) is one of the top 4 provinces in the country, which explains why so many listings are in this province.\n",
    "\n",
    "- **Bs.As. G.B.A. Zona Norte (Northern Zone of Greater Buenos Aires)**\n",
    "  - **Class Type**: Affluent, upper-middle to upper-class.\n",
    "  - **Real Estate**: High property values, upscale residential areas, and gated communities.\n",
    "\n",
    "- **Capital Federal (Autonomous City of Buenos Aires)**\n",
    "  - **Class Type**: Diverse, ranging from working-class to upper-class.\n",
    "  - **Real Estate**: Diverse property market with high values in central neighborhoods like Palermo, Recoleta, and Puerto Madero.\n",
    "\n",
    "- **Bs.As. G.B.A. Zona Oeste (Western Zone of Greater Buenos Aires)**\n",
    "  - **Class Type**: Mixed, predominantly middle-income and working-class.\n",
    "  - **Real Estate**: Generally lower property values compared to the northern zone, with a mix of residential and industrial areas.\n",
    "\n",
    "- **Bs.As. G.B.A. Zona Sur (Southern Zone of Greater Buenos Aires)**\n",
    "  - **Class Type**: Predominantly middle to lower-income.\n",
    "  - **Real Estate**: Lower property values, with some areas undergoing urban renewal and modernization projects.\n",
    "\n",
    "#### **Investigate lat-lon**\n",
    "\n",
    "Next, we will extract the latitude and longitude coordinates from the `lat-lon` column. \n",
    "\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2165e03f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split the `lat-lon` column into multiple columns\t\n",
    "df_buenos_aires[['latitude', 'longitude']] = df_buenos_aires['lat-lon'].str.split(',', expand=True).astype(float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f13e9525",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.scatter_mapbox( \n",
    "    df_buenos_aires.dropna(subset=['latitude', 'longitude']),  # Our DataFrame without missing values in `latitude` and `longitude`\n",
    "    lat='latitude',\n",
    "    lon='longitude',\n",
    "    width=900,  # Width of map\n",
    "    height=800,  # Height of map\n",
    "    color='price_aprox_usd',\n",
    "    hover_data=[\"price_aprox_usd\"],  # Display price when hovering mouse over house\n",
    ")\n",
    "\n",
    "fig.update_layout(mapbox_style=\"open-street-map\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d9a51e6",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "From the map, we can see that most listings are located in the `Capital Federal` province of Buenos Aires, specifically in the`Autonomous City of Buenos Aires`.\n",
    "\n",
    "We also observe that we have some listings that are not in Argentina. Possibly, it was due to the scraping process.\n",
    "We will focus on the listings in the `Capital Federal` province. Below is the map of the `Capital Federal` province and all the departments in the city.\n",
    "\n",
    "<img src= \"https://www.mapsland.com/maps/south-america/argentina/buenos-aires/large-tourist-and-administrative-subdivisions-map-of-buenos-aires-city.jpg\" width=60% style=\"display: block; margin-left: auto; margin-right: auto;\"/>\n",
    "\n",
    "\n",
    "## **Data preparation (preprocessing)**\n",
    "\n",
    "### **Data cleaning**\n",
    "\n",
    "Now, we will proceed with the data cleaning process and we will base our process on the EDA results.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be68377a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Columns with more than 50% of values missing\t\n",
    "def nan_columns(df):\n",
    "    \"\"\"Columns with more than 50% of values missing\"\"\"\n",
    "    cols_series = df.isna().sum() / len(df)\n",
    "    return (cols_series[cols_series >= 0.5]).index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1616df5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def wrangle(filepath):\n",
    "    \"\"\"\n",
    "    Read and wrangle data from CSV file\n",
    "\n",
    "    Parameters\n",
    "    ---------- \n",
    "    filepath: str\n",
    "        Path to CSV file\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    df: pandas.DataFrame\n",
    "        Wrangled data\n",
    "    \"\"\"\n",
    "    \n",
    "    # Read CSV file\n",
    "    df = pd.read_csv(filepath, low_memory=False)\n",
    "\n",
    "    # Subset data: Apartments in \"Capital Federal\", less than $400,000 \n",
    "    mask_cf = df[\"place_with_parent_names\"].str.contains(\"capital federal\", case=False)\n",
    "    mask_apt = df[\"property_type\"] == \"apartment\"\n",
    "    mask_price = df[\"price_aprox_usd\"] <= 400_000\n",
    "    df = df[mask_cf & mask_apt & mask_price]\n",
    "\n",
    "    # Subset data: Remove outliers for \"surface_total_in_m2\" and \"surface_covered_in_m2\"\n",
    "    # Values are based on the get_quantile_category() analysis\n",
    "    mask_sc = df['surface_covered_in_m2'].between(34, 500)\n",
    "    mask_st = df['surface_total_in_m2'] > 0\n",
    "    df = df[mask_sc & mask_st]\n",
    "\n",
    "    # Split \"lat-lon\" column\n",
    "    df[[\"latitude\", \"longitude\"]] = df[\"lat-lon\"].str.split(\",\", expand=True).astype(float)\n",
    "\n",
    "    # Split \"place_with_parent_names\" column\n",
    "    df['municipality'] = df.place_with_parent_names.str.split('|', expand=True)[3]\n",
    "    \n",
    "    # Drop unnecessary columns\n",
    "    df.drop(columns=[\"place_with_parent_names\", \"lat-lon\", \"price_per_m2\",], axis=1, inplace=True)\n",
    "    \n",
    "    # Drop null columns (50%+ null rows)\n",
    "    ## Namely `floor`, `expenses` and `price_usd_per_m2`\n",
    "    nan_cols = nan_columns(df)\n",
    "    df.drop(columns=nan_cols, inplace=True)\n",
    "\n",
    "    # Drop low- and high-cardinality categorical variables\n",
    "    df.drop(columns=['operation', 'property_type', 'currency', 'properati_url'], inplace=True)\n",
    "\n",
    "    # Drop leaky columns\n",
    "    df.drop(columns=['price', 'price_aprox_local_currency'], inplace=True)\n",
    "\n",
    "    # Remove empty rows from municipality\n",
    "    df = df[df['municipality'].str.len() > 0]\n",
    "\n",
    "    # Proper headings for our data\n",
    "    proper_headings = ['municipality', 'latitude', 'longitude', 'surface_total_in_m2', 'surface_covered_in_m2', 'rooms', 'price_aprox_usd']\n",
    "    \n",
    "    return df[proper_headings].reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b8b539f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run wrangle function\t\n",
    "df = wrangle('data/buenos-aires-real-estate.csv')\n",
    "\n",
    "# Create a copy of the original dataframe for later use\n",
    "df_full = df.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e729d5c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc7668c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb89fa39",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "Our data is clean and ready to use. We filtered the data by province to only focus on the `Capital Federal` province. We also dropped the columns with more than 50% of values missing and set the limit for the price to be less than $400,000. We also split the `lat-lon` column into two columns `latitude` and `longitude` and dropped the original column. \n",
    "\n",
    "We also dropped the columns `price_aprox_local_currency`, `currency`, `price_aprox_usd`, and `price` because we observed inconsistent data. We also dropped the columns `price_per_m2` and `price_usd_per_m2` because they were leaky. \n",
    "\n",
    "\n",
    "### **Visualize Data**\n",
    "\n",
    "We will visualize the data after wrangling to gain insights. We will start by examining the average price of the houses by municipality.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9301fcf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Horizontal bar chart of Average price by municipality\n",
    "plt.figure(figsize=(10, 13))\n",
    "(df.groupby('municipality')\n",
    " .price_aprox_usd.mean().sort_values(ascending=True)\n",
    " .plot.barh(color='royalblue', alpha=0.8, edgecolor='black')\n",
    " )\n",
    "\n",
    "plt.ylabel('')\n",
    "plt.xlabel('Average price in USD')\n",
    "plt.title('Average Apartment price by municipality\\n')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6b8df4d",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "From the above chart, we conclude the following:\n",
    "\n",
    "- **High Prices in the North and East:**\n",
    "    - The northern and eastern neighborhoods of Buenos Aires, especially those close to the city center and\n",
    "    waterfront, tend to have the highest apartment prices. This is likely due to their proximity to business districts, cultural attractions, and better infrastructure.\n",
    "\n",
    "- **Moderate Prices in Central Areas:**\n",
    "    - Central neighborhoods have moderate to high prices, reflecting a mix of residential and commercial use, with good access to amenities and transport.\n",
    "\n",
    "- **Lower Prices in the South and West:**\n",
    "    - The southern and western parts of Buenos Aires generally have lower apartment prices. These areas are more residential, with some industrial zones and lower socio-economic conditions compared to the north and east.\n",
    "\n",
    "\n",
    "Moving on, we will explore the general trend of price vs. total surface and price vs. covered surface.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "810b1f18",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scatter plot of price vs. total_surface_in_m2 and price vs. covered_suface_in_m2 one next to the other \n",
    "# You can zoom in and out \n",
    "\n",
    "# Price vs. total surface\n",
    "fig = px.scatter(\n",
    "    df,\n",
    "    x=\"surface_total_in_m2\",\n",
    "    y=\"price_aprox_usd\",\n",
    "    color=\"price_aprox_usd\",\n",
    "    title=\"Price vs. total surface in m\\u00B2\",\n",
    "    labels={\"surface_total_in_m2\": \"Total surface in m\\u00B2\", \"price_aprox_usd\": \"Price in USD\"},\n",
    "    template=\"plotly_white\",\n",
    "    width=800,\n",
    "    height=600\n",
    ")\n",
    "\n",
    "fig.update_traces(\n",
    "    marker=dict(size=8, line=dict(width=1, color='white')),\n",
    "    selector=dict(mode='markers')\n",
    "    )\n",
    "\n",
    "fig.show()\n",
    "\n",
    "# Price vs. covered surface\n",
    "fig = px.scatter(\n",
    "    df,\n",
    "    x=\"surface_covered_in_m2\",\n",
    "    y=\"price_aprox_usd\",\n",
    "    color=\"price_aprox_usd\",\n",
    "    title=\"Price vs. covered surface in m\\u00B2\",\n",
    "    labels={\"suface_covered_in_m2\": \"Covered surface in m\\u00B2\", \"price_aprox_usd\": \"Price in USD\"},\n",
    "    template=\"plotly_white\",\n",
    "    width=800,\n",
    "    height=600\n",
    ")\n",
    "\n",
    "fig.update_traces(\n",
    "    marker=dict(size=8, line=dict(width=1, color='white')),\n",
    "    selector=dict(mode='markers')\n",
    "    )\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d37c785",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "We can see that there is a positive correlation between price and surface area. When surface area is higher, the price tends to increase. This holds true for both `total_surface_in_m2` and `covered_suface_in_m2`.\n",
    "\n",
    "There is a slight curvature toward the top of the scatter plot in the case of `total_surface_in_m2` which suggests an exponential relationship. Inversely, there is a slight curvature toward the bottom of the plot in the case of `covered_suface_in_m2`, which suggests a logarithmic relationship.\n",
    "\n",
    "Examining the spread of the data, we can see that variability is more pronounced in the `covered_surface_in_m2` column than in the `total_suface_in_m2` column.\n",
    "\n",
    "It might not be clearly visible, but with the help of zoom and box selection, we can observe the phenomenon.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e05b69eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot 3d scatter plot of price vs. latitude and longitude \n",
    "# With samll markers and color gradient\n",
    "\n",
    "fig = px.scatter_3d(\n",
    "    df,\n",
    "    x=\"latitude\",\n",
    "    y=\"longitude\",\n",
    "    z=\"price_aprox_usd\",\n",
    "    color=\"price_aprox_usd\",\n",
    "    title=\"Price vs. latitude and longitude\",\n",
    "    labels={\"longitude\": \"Longitude\", \"latitude\": \"Latitude\", \"price_aprox_usd\": \"Price in USD\"},\n",
    "    template=\"plotly_white\",\n",
    "    width=800,\n",
    "    height=600\n",
    ")   \n",
    "\n",
    "fig.update_traces(\n",
    "    marker=dict(size=4, line=dict(width=1, color='white')),\n",
    "    selector=dict(mode='markers')\n",
    "    )\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f429ad1",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "From the 3D scatter plot, we can see that as we move towards the coast, the price tends to increase.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38e561c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scatter plot of price vs. latitude and longitude with price as a color gradient\n",
    "# Map is open-street-map and zoomed in\n",
    "fig = px.scatter_mapbox( \n",
    "    df,\n",
    "    lat='latitude',\n",
    "    lon='longitude',\n",
    "\n",
    "    mapbox_style=\"open-street-map\",  # Mapbox style\n",
    "\n",
    "    zoom=11,  # Zoom level\n",
    "    title=\"Price vs. latitude and longitude\",\n",
    "    labels={\"longitude\": \"Longitude\", \"latitude\": \"Latitude\", \"price_aprox_usd\": \"Price in USD\"},\n",
    "\n",
    "    width=900,  # Width of map\n",
    "    height=800,  # Height of map\n",
    "\n",
    "    color='price_aprox_usd',\n",
    "    size='price_aprox_usd',\n",
    "    \n",
    "    hover_data=[\"price_aprox_usd\"],  # Display price when hovering mouse over house\n",
    ")\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8b29528",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "The Scatter map plot is showing clearly what the 3D Scatter plot is showing. As we move towards the coast, the price tends to increase.\n",
    "\n",
    "## **Modeling**\n",
    "\n",
    "We proceed by specifying features and target variable and then Split the data into training and testing sets. We hold 20% for the testing set.\n",
    "\n",
    "We notice that `municipality` is an **out-of-sample feature**. This happens when certain features (like municipalities) present in the test set were not seen during training, the model may not have learned the necessary patterns to make accurate predictions for these unseen categories. \n",
    "\n",
    "Our decision is the following, instead of relying solely on municipality names, we will incorporate latitude and longitude directly as features. This allows our model to learn spatial relationships which might help in generalizing to unseen municipalities.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "475bacfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop the `department_municipality` column\n",
    "df.drop(columns=['municipality'], inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c753aa0",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "We will examine the correlation between features and target variable after cleaning the data and after dropping the `municipality` column.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0092cae6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Correlation matrix between numerical columns\n",
    "corr_frame = df.corr()  # type: ignore\t\n",
    "\n",
    "# Mask the upper triangle of the correlation matrix\n",
    "mask_corr_df = np.triu(np.ones_like(corr_frame, dtype=bool))\n",
    "\n",
    "plt.figure(figsize = (12, 8))\n",
    "sns.heatmap(corr_frame, annot = True, cmap = \"Spectral\", linewidths = 2, linecolor = \"white\", fmt='.2f', mask=mask_corr_df) # type: ignore\n",
    "plt.title(\"Correlation Matrix\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8feec3df",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "First, we notice that `total_surface_in_m2` and `covered_suface_in_m2` are strongly correlated. This can bring high multi-collinearity into the model. We will have to drop one of them. Our choice is the `total_surface_in_m2` feature.\n",
    "\n",
    "Next, we notice that `price_aprox_usd` is highly correlated with `total_surface_in_m2` and `covered_suface_in_m2`. Also, `rooms` is  moderately correlated with `price_aprox_usd` and highly correlated with `total_surface_in_m2` and `covered_suface_in_m2`. \n",
    "\n",
    "In summary, we will drop `total_surface_in_m2` and `rooms` features.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71937f8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop `rooms` and `surface_total_in_m2` columns\n",
    "\n",
    "df.drop(columns=['rooms', 'surface_total_in_m2'], axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77718fb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "corr_frame_drop = df.corr()  # type: ignore\n",
    "\n",
    "# Mask the upper triangle of the correlation matrix\n",
    "mask_corr_df_drop = np.triu(np.ones_like(corr_frame_drop, dtype=bool))\n",
    "\n",
    "plt.figure(figsize = (10, 8))\n",
    "sns.heatmap(corr_frame_drop, annot = True, cmap = \"Spectral\", linewidths = 2, linecolor = \"white\", fmt='.2f', mask=mask_corr_df_drop) # type: ignore\n",
    "plt.title(\"Correlation Matrix \\n After Dropping `rooms` and `total_surface_in_m2`\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5678cf4a",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "### **Split the data**\n",
    "\n",
    "We will split the data into training and testing sets with 20% of data in the test set, and then scale the data to have zero mean and unit standard deviation.\n",
    "First, we will check if there are some missing values in the data.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5e71b98",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.isna().sum(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "091f84a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.dropna(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cec79469",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"We have {df.shape[0]:,} rows and {df.shape[1]} columns.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2c99d1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify features and target\n",
    "\n",
    "X = df.drop(columns=['price_aprox_usd'], axis=1)\n",
    "y = df.price_aprox_usd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f66bcac",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split data into train and test sets with 20% of data in test set\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d702049c",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('After split:')\n",
    "print(f\"We have {X_train.shape[0]:,} rows in the training set and {X_test.shape[0]:,} rows in the test set.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5681244",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Baseline model\n",
    "y_pred_baseline = np.full(y_train.shape[0], y_train.mean())\n",
    "\n",
    "mae_baseline = mean_absolute_error(y_train, y_pred_baseline) \n",
    "\n",
    "print(f\"Baseline MAE: ${mae_baseline:,.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b2ba3d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "(59897.91 - 23809.9879) * 100 / 59897.91"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93c602d1",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "We'll have to beat the baseline model by at least 10% in terms of mean absolute error.\n",
    "\n",
    "### **Iterate**\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f996584f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add constant to X_train, X_test \n",
    "X_train = sm.add_constant(X_train)  \n",
    "X_test = sm.add_constant(X_test)\n",
    "\n",
    "# Fit the model\n",
    "model_0 = sm.OLS(y_train, X_train)\n",
    "fitted_0 = model_0.fit()\n",
    "\n",
    "# Print summary\t\n",
    "print(fitted_0.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8f03256",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred_ols_0 = fitted_0.predict(X_test) \n",
    "\n",
    "# MAE before inverse scaling (bis)\n",
    "mae_ols_0 = mean_absolute_error(y_true=y_test, y_pred=y_pred_ols_0)\n",
    "print(f\"OLS MAE: ${mae_ols_0:,.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9839b0da",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare MAE\n",
    "pct_decrease_mae_0 = (mae_baseline - mae_ols_0) * 100 / mae_baseline \n",
    "print(f\"OLS MAE Decrease: {pct_decrease_mae_0:.2f}%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "852a18b0",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "Great! We beat the baseline model by 45.6% in terms of mean absolute error. \n",
    "\n",
    "#### **Detect Heteroskedasticity using formal test**\n",
    "\n",
    "\n",
    "The null hypothesis $( H_0 )$ in the context of heteroscedasticity tests is typically formulated as follows:\n",
    "\n",
    "$H_0: \\text{The errors/residuals exhibit homoscedasticity (constant variance).}$\n",
    "\n",
    "In simpler terms, the null hypothesis states that there is no systematic difference in the variance of errors across different levels or values of the independent variables. If the p-value from the heteroscedasticity test is sufficiently small, we reject this null hypothesis in favor of the alternative hypothesis, indicating that there is evidence of heteroscedasticity present in the data.\n",
    "To determine whether to reject or not reject the null hypothesis $( H_0)$ based on the p-value:\n",
    "\n",
    "1. **Reject $( H_0 )$ (Evidence of Effect):**\n",
    "   - Small p-value (p < α, where α is the significance level): Strong evidence against $( H_0 )$.\n",
    "   - Indicates that observed data is unlikely under the assumption of $( H_0 )$.\n",
    "\n",
    "2. **Do Not Reject $( H_0 )$ (Lack of Evidence of Effect):**\n",
    "   - Large p-value (p ≥ α): Not enough evidence to conclude against $( H_0 )$.\n",
    "   - Suggests that observed data is consistent with $( H_0 )$ or does not deviate significantly from expected values under $( H_0 )$.\n",
    "\n",
    "In the context of heteroscedasticity tests:\n",
    "- Small p-value: Evidence of heteroscedasticity; reject the assumption of constant variance.\n",
    "- Large p-value: Insufficient evidence of heteroscedasticity; do not reject the assumption of constant variance.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a264fca4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.stats.diagnostic import het_breuschpagan\n",
    "\n",
    "# Perform Breusch-Pagan test\n",
    "bp_test = het_breuschpagan(fitted_0.resid, fitted_0.model.exog) # type: ignore\n",
    "print(\"Breusch-Pagan test p-value:\", bp_test[1])\n",
    "\n",
    "# Check if p-value is less than the significance level\n",
    "alpha = 0.05\n",
    "if bp_test[1] < alpha:\n",
    "    print(f\"Reject the null hypothesis: Heteroskedasticity is present.\")\n",
    "else:\n",
    "    print(f\"Fail to reject the null hypothesis: Heteroskedasticity is not present.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e085e94d",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "We got a p-value of **7.74e-80**, which exteremely small. Therefore, we reject the null hypothesis $( H_0 )$ in favor of the alternative hypothesis, indicating that there is evidence of heteroscedasticity present in the data.\n",
    "\n",
    "#### **Detect Heteroskedasticity using visual methods: Residuals plot**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6b0c243",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot residuals vs. predicted values\n",
    "\n",
    "plt.figure(figsize=(10, 8))\n",
    "sns.scatterplot(x=y_train, y=fitted_0.resid, color='royalblue', edgecolor='black', alpha=0.8) \n",
    "plt.axhline(y=0, color='r', linestyle='--', linewidth=1)\n",
    "plt.title(\"Residuals vs. Predicted Values\\n (On Training Set)\")\n",
    "plt.xlabel(\"Predicted Values\")\n",
    "plt.xticks(rotation=0)\n",
    "plt.ylabel(\"Residuals\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0f89db2",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif;\">\n",
    "\n",
    "We can see that as predicted values increase, the variability of the residuals increases too, indicating heteroscedasticity.\n",
    "\n",
    "We will deal with it by transforming the target variable and some of the features as needed..\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7820ccae",
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.graphics.regressionplots import plot_leverage_resid2\n",
    "\n",
    "# plot leverage vs. residuals\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 8))\n",
    "plot_leverage_resid2(fitted_0, ax=ax, markersize=10, markeredgecolor=\"black\", markerfacecolor=\"royalblue\")  # type: ignore\n",
    "# Customize text properties\n",
    "for text in ax.texts:\n",
    "    # uncomment the following line if you want to adjust the text size\n",
    "    # text.set_fontsize(8)\n",
    "\n",
    "    # set visibility to True if you want to show the observations indices on the plot\n",
    "    text.set_visible(False)  \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14d90124",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "The graph above shows the leverage vs. normalized residuals. We can see that the leverage is high for some of the observations, which may lead to our model not being able to predict accurately. High leverage points are considered influential in the model.\n",
    "Testing the assumption of **Normality of residuals**:\n",
    "\n",
    "The residuals should exhibit a normal distribution. This is easiest to detect by inspecting a QQ-plot of residuals and a residuals vs dependent variable plot.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2e548ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot sample quantiles vs. the theoretical quantiles \n",
    " \n",
    "fig = sm.graphics.qqplot(fitted_0.resid, dist=sps.norm, line='45', fit=True, marker=\"o\", markersize=8, markeredgecolor=\"black\", markerfacecolor=\"royalblue\")\n",
    "plt.title(\"Sample Quantiles vs. Theoretical Quantiles\", fontsize=18, fontweight=\"bold\") \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33cee969",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "We can see that the residuals are not normally distributed. The second bisection of the QQ-plot indicates the case where the residuals are normally distributed, and from examining the graph we spot some a trend not coresponding to the normal distribution. \n",
    "\n",
    "## **Reverse Heteroskedasticity**\n",
    "\n",
    "### **Detect Influential Points**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2996b09",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the influence statistics\n",
    "influence = fitted_0.get_influence() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea16e2a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Transform the influence statistics into a DataFrame\n",
    "influence_df = influence.summary_frame()\n",
    "influence_df.head(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3ae1e55",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "We will rely on the `Summary Frame` to guide us in our analysis in reversing heteroskedasticity.\n",
    "\n",
    "#### **Visual Ispection**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c747271d",
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold_leverage = 3 * influence_df.cooks_d.mean()\n",
    "number_influential = influence_df[influence_df.cooks_d > threshold_leverage].shape[0]\n",
    "\n",
    "print(f\"Number of Influential points with the threshold of {threshold_leverage}, is {number_influential}.\")\n",
    "\n",
    "plt.figure(figsize=(10, 8))\n",
    "plt.stem(np.arange(len(influence_df)), influence_df.cooks_d, markerfmt=\",\")\n",
    "plt.axhline(y=threshold_leverage, color='r', linestyle='--')\n",
    "plt.ylim(0, 0.05)\n",
    "plt.xlabel('Observation Index')\n",
    "plt.ylabel('Cook\\'s Distance')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0de7995",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "`Cook's Distance` is a measure of the influence of an observation on the regression model. We can see that there are more than 15 influential points with a Cook's Distance above the threshold chosen. We will have to find the sweet spot to remove these influential points and at the same time reverse the heteroskedasticity effect.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b7bf6d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold_residual = 2.5\n",
    "\n",
    "plt.figure(figsize=(10, 8))\n",
    "plt.stem(np.arange(len(influence_df)), influence_df.standard_resid, markerfmt=\",\")\n",
    "plt.axhline(y=0, color='r', linestyle='--')\n",
    "plt.axhline(y=-threshold_residual, color='r')\n",
    "plt.axhline(y=threshold_residual, color='r')\n",
    "\n",
    "plt.xlabel('Observation Index')\n",
    "plt.ylabel('Standard Residuals')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc2ab00c",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "The `Standard Residuals` vs. `Observation Index` plot shows that there are some outliers in the data based on the threshold chosen.\n",
    "\n",
    "#### **Remove high leverage & high residual points**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c45e7695",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "high_leverage_threshold = threshold_leverage\n",
    "high_residual_threshold = threshold_residual\n",
    "\n",
    "# Identify high leverage and high residual points\n",
    "high_leverage_points = influence_df.cooks_d > high_leverage_threshold\n",
    "high_residual_points = np.abs(influence_df.standard_resid) > high_residual_threshold\n",
    "\n",
    "# Combine both conditions to find points with either high leverage or high residuals\n",
    "influential_points = high_leverage_points & high_residual_points\n",
    "\n",
    "# Print the indices of the influential points\n",
    "print(f\"We found {len(np.where(influential_points)[0])} influential points that has high leverage and high residuals.\")\n",
    "print(\"Indices of influential points:\\n\", np.where(influential_points)[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3489d7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove influential points from the Train set\n",
    "filtered_X_train = X_train[~influential_points]\n",
    "filtered_y_train = y_train[~influential_points]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f956342",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "#### **Fit model on Filtered Data**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "687299c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit the model on the filtered Train set and print summary\n",
    "model_1 = sm.OLS(filtered_y_train, filtered_X_train)\n",
    "fitted_1 = model_1.fit()\n",
    "\n",
    "print(fitted_1.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e50bbaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot sample quantiles vs. the theoretical quantiles after transformation\n",
    " \n",
    "fig = sm.graphics.qqplot(fitted_1.resid, dist=sps.norm, line='45', fit=True, marker=\"o\",\n",
    "                          markersize=8, markeredgecolor=\"black\", markerfacecolor=\"royalblue\")\n",
    "plt.title(\"Sample Quantiles vs. Theoretical Quantiles\", fontsize=18, fontweight=\"bold\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b227e7ab",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "The Q-Q plot shows that the residuals are slightly now normally distributed. With some inconsistent behaviour in upper quantiles.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b359acc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Predict on the test set\n",
    "y_pred_ols_1 = fitted_1.predict(X_test)\n",
    "\n",
    "# Calculate MAE for the test set\n",
    "mae_ols_1 = mean_absolute_error(y_true=y_test, y_pred=y_pred_ols_1)\n",
    "print(f\"OLS MAE: ${mae_ols_1:,.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76a5dd91",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "The MAE has decreased slightly when we removed influential points.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8910d982",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare MAE & RMSE\n",
    "\n",
    "pct_decrease_mae_1 = (mae_baseline - mae_ols_1) * 100 / mae_baseline \n",
    "\n",
    "print(f\"OLS MAE Decrease (After): {pct_decrease_mae_1:.2f}%\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35193a42",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform Breusch-Pagan test\n",
    "bp_test = het_breuschpagan(fitted_1.resid, fitted_1.model.exog) # type: ignore\n",
    "print(\"Breusch-Pagan test p-value:\", bp_test[1])\n",
    "\n",
    "# Check if p-value is less than the significance level\n",
    "alpha = 0.05\n",
    "if bp_test[1] < alpha:\n",
    "    print(f\"Reject the null hypothesis: Heteroskedasticity is present.\")\n",
    "else:\n",
    "    print(f\"Fail to reject the null hypothesis: Heteroskedasticity is not present.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d7e847c",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "Although the model has improved slightly, the heteroskedasticity is still present. This suggest that we need some advanced techniques to deal with it.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33599dd8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot residuals vs. predicted values\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "sns.scatterplot(x=y_train, y=fitted_1.resid, color='royalblue', edgecolor='black', alpha=0.8) \n",
    "plt.axhline(y=0, color='r', linestyle='--', linewidth=1)\n",
    "plt.title(\"Residuals vs. Predicted Values\\n (After)\")\n",
    "plt.xlabel(\"Predicted Values\")\n",
    "plt.xticks(rotation=0)\n",
    "plt.ylabel(\"Residuals\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9c92669",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "The range for the `Predicted Values` has chrunked due to influential point removal. But variance of the residuals is still not constant.\n",
    "Some of the advanced techniques to deal with heteroskedasticity are **Box-Cox transformation**. We will use this technique to transform the both the target variable and the `surface_covered_in_m2` variable, because from previous analysis, we saw that `surface_covered_in_m2` is right skewed. \n",
    "\n",
    "#### **Transform Data: Box-Cox**\n",
    "We will build a function to apply Box-Cox transformation to both the target variable and the `surface_covered_in_m2` variable and choose the best lambda value that will minimize the MAE error and invert the heteroskedasticity effect.\n",
    "\n",
    "We will rely on the `Breusch-Pagan test` to determine if heteroskedasticity is present.\n",
    "##### **Find best lambda**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85f9f0be",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import  StandardScaler\n",
    "from scipy.special import inv_boxcox\n",
    "\n",
    "# Define lambda values to test\n",
    "def test_lambda(lambdas, filtered_X_train, filtered_y_train, X_test, y_test):\n",
    "\n",
    "    bp_list = list()\n",
    "    mae_list = list()\n",
    "\n",
    "    for lambda_ in lambdas:\n",
    "        X_train_tl = filtered_X_train.copy()\n",
    "        X_test_tl = X_test.copy()\n",
    "\n",
    "        y_train_tl = filtered_y_train.copy()\n",
    "        y_test_tl = y_test.copy()\n",
    "        \n",
    "        # Apply Box-Cox transformation to `covered_suface_in_m2`\n",
    "        X_train_tl.surface_covered_in_m2 = sps.boxcox(X_train_tl.surface_covered_in_m2, lambda_)\n",
    "        X_test_tl.surface_covered_in_m2 = sps.boxcox(X_test_tl.surface_covered_in_m2, lambda_)\n",
    "\n",
    "        # Apply Box-Cox transformation to `price_usd`\n",
    "        y_train_tl = sps.boxcox(y_train_tl, lambda_)\n",
    "        y_test_tl = sps.boxcox(y_test_tl, lambda_)\n",
    "\n",
    "        # Fit the model\n",
    "        model = sm.OLS(y_train_tl, X_train_tl) \n",
    "        fitted = model.fit()\n",
    "        \n",
    "        # Perform Breusch-Pagan test\n",
    "        bp_test = het_breuschpagan(fitted.resid, fitted.model.exog)\n",
    "        bp_list.append(bp_test[1])\n",
    "\n",
    "        # Predict on the test set\n",
    "        y_pred_ols_test = fitted.predict(X_test_tl)\n",
    "\n",
    "        # Invert the Box-Cox transformation by s.inv_boxcox\n",
    "        y_pred_ols_test_inv = inv_boxcox(y_pred_ols_test, lambda_)\n",
    "\n",
    "        # Calculate MAE\n",
    "        mae_ols = mean_absolute_error(y_true=y_test, y_pred=y_pred_ols_test_inv)\n",
    "        mae_list.append(mae_ols)\n",
    "\n",
    "    df_results = pd. DataFrame({'lambda_value': lambdas, 'Breusch-Pagan test p-value': bp_list, 'MAE': mae_list})\n",
    "\n",
    "    # Filter df_results for p-value > 0.05 which means we fail to reject the null hypothesis that Heteroskedasticity is present\n",
    "    df_results = df_results[df_results['Breusch-Pagan test p-value'] > 0.05]  # Heteroskedasticity is not present \n",
    "\n",
    "    return df_results.sort_values(by='MAE', ascending=True).reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1d57cc9",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "We will apply `test_lambda` function to `X_train`, `y_train`, `X_test`, `y_test` and feed the function with `np.arange(-0.5, 0.5, 0.01)` as the range for the possible lambda values we believe will minimize the MAE error.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99dd17de",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_lambda_df = test_lambda(np.arange(-0.5, 0.5, 0.01), X_train, y_train, X_test, y_test)\n",
    "test_lambda_df[:5] # top 5 rows"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6c6d446",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "The results show that the best lambda value is **-0.03**, and as we can see the p-value is greater than 0.05, so we fail to reject the null hypothesis that Heteroskedasticity is present.\n",
    "\n",
    "##### **Fit model**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "831d08be",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_tr = X_train.copy()\n",
    "y_train_tr = y_train.copy()\n",
    "\n",
    "# Retrieve the best lambda value from the test_lambda_df\n",
    "best_lambda = test_lambda_df.lambda_value[0]\n",
    "\n",
    "# Apply Box-Cox transformation to `covered_suface_in_m2` & `price_aprox_usd` using the best lambda value\n",
    "X_train_tr.surface_covered_in_m2 = sps.boxcox(X_train_tr.surface_covered_in_m2, best_lambda) # type: ignore\n",
    "y_train_tr = sps.boxcox(y_train_tr, best_lambda)\n",
    "\n",
    "# Fit the model on the filtered Train set and print summary\n",
    "model_2 = sm.OLS(y_train_tr, X_train_tr)\n",
    "fitted_2 = model_2.fit()\n",
    "\n",
    "print(fitted_2.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76edee59",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform Breusch-Pagan test\n",
    "bp_test = het_breuschpagan(fitted_2.resid, fitted_2.model.exog) # type: ignore\n",
    "print(\"Breusch-Pagan test p-value:\", bp_test[1])\n",
    "\n",
    "# Check if p-value is less than the significance level\n",
    "alpha = 0.05\n",
    "if bp_test[1] < alpha:\n",
    "    print(f\"Reject the null hypothesis: Heteroskedasticity is present.\")\n",
    "else:\n",
    "    print(f\"Fail to reject the null hypothesis: Heteroskedasticity is not present.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53cb5ea9",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_test_tr = X_test.copy()\n",
    "\n",
    "X_test_tr.surface_covered_in_m2 = sps.boxcox(X_test_tr.surface_covered_in_m2, best_lambda)  # type: ignore\n",
    "\n",
    "# Predict on the transformed test set\n",
    "y_pred_ols_2 = fitted_2.predict(X_test_tr)\n",
    "\n",
    "\n",
    "# Invert the Box-Cox transformation by inv_boxcox\n",
    "y_pred_ols_2_inv = inv_boxcox(y_pred_ols_2, best_lambda)\n",
    "\n",
    "\n",
    "# Calculate MAE for the test set\n",
    "mae_ols_2 = mean_absolute_error(y_true=y_test, y_pred=y_pred_ols_2_inv)\n",
    "print(f\"OLS MAE: ${mae_ols_2:,.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96554dd7",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "As we anticipated, the MAE has decreased from iteration 1 to iteration 2. \n",
    "\n",
    "##### **Transform Filtered Data: Box-Cox**\n",
    "\n",
    "Now, we will apply Box-Cox transformation to `price_aprox_usd` & `covered_surface_in_m2` on the filtered train and test sets.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e81edd0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply Box-Cox transformation to `price_usd` & `covered_surface_in_m2`\n",
    "X_train_sc = filtered_X_train.copy()\n",
    "y_train_sc = filtered_y_train.copy()\n",
    "\n",
    "X_test_sc = X_test.copy()\n",
    "y_test_sc = y_test.copy()\n",
    "\n",
    "best_lambda = test_lambda_df.lambda_value[0]\n",
    "X_train_sc.surface_covered_in_m2 = sps.boxcox(X_train_sc.surface_covered_in_m2, best_lambda)  # type: ignore\n",
    "y_train_sc = sps.boxcox(y_train_sc, best_lambda)\n",
    "\n",
    "X_test_sc.surface_covered_in_m2 = sps.boxcox(X_test_sc.surface_covered_in_m2, best_lambda)  # type: ignore\n",
    "y_test_sc = sps.boxcox(y_test_sc, best_lambda)\n",
    "\n",
    "# Scale Data\n",
    "std_scaler = StandardScaler(with_mean=True, with_std=True, copy=True)\n",
    "\n",
    "X_train_sc = pd.DataFrame(std_scaler.fit_transform(X_train_sc), columns=X_train.columns, index=filtered_X_train.index)  # type: ignore\n",
    "X_test_sc = pd.DataFrame(std_scaler.transform(X_test_sc), columns=X_test.columns, index=X_test.index)  # type: ignore\n",
    "X_train_sc.const = 1\n",
    "X_test_sc.const = 1\n",
    "\n",
    "\n",
    "# Fit the model on the filtered Train set and print summary\n",
    "model_3 = sm.OLS(y_train_sc, X_train_sc)\n",
    "fitted_3 = model_3.fit()\n",
    "\n",
    "print(fitted_3.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8e1d334",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform Breusch-Pagan test\n",
    "bp_test = het_breuschpagan(fitted_3.resid, fitted_3.model.exog) # type: ignore\n",
    "print(\"Breusch-Pagan test p-value:\", bp_test[1])\n",
    "\n",
    "# Check if p-value is less than the significance level\n",
    "alpha = 0.05\n",
    "if bp_test[1] < alpha:\n",
    "    print(f\"Reject the null hypothesis: Heteroskedasticity is present.\")\n",
    "else:\n",
    "    print(f\"Fail to reject the null hypothesis: Heteroskedasticity is not present.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec33269d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Predict on the transformed test set\n",
    "y_pred_ols_3 = fitted_3.predict(X_test_sc)\n",
    "\n",
    "# Invert the Box-Cox transformation by inv_boxcox\n",
    "y_pred_ols_3_inv = inv_boxcox(y_pred_ols_3, best_lambda)\n",
    "\n",
    "# Calculate MAE for the test set\n",
    "mae_ols_3 = mean_absolute_error(y_true=y_test, y_pred=y_pred_ols_3_inv)\n",
    "print(f\"OLS MAE: ${mae_ols_3:,.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c72dd91",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "Although the model R-squared and MAE is improved, the heteroskedasticity is still present. Our model has a very tight interval for `leverage` and `residuals` thresholds to exclude influential points and reverse heteroskedasticity.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c59346f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot residuals vs. predicted values\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "sns.scatterplot(x=y_train, y=fitted_3.resid, color='royalblue', edgecolor='black', alpha=0.8) \n",
    "plt.axhline(y=0, color='r', linestyle='--', linewidth=1)\n",
    "plt.title(\"Residuals vs. Predicted Values\\n (After)\")\n",
    "plt.xlabel(\"Predicted Values\")\n",
    "plt.xticks(rotation=0)\n",
    "plt.ylabel(\"Residuals\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6969ec6",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "We still observe heteroskedasticity in the residuals.\n",
    "\n",
    "We will use another algorithm on our filtered dataset, which will incorporate all the insights we have gathered so far from the OLS model.\n",
    "\n",
    "We have chosen `Gradient Boosting Regressor` as the final model.\n",
    "\n",
    "## **Final Model: Gradient Boosting Regressor**\n",
    "\n",
    "### **Fiting the model on cleaned data**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dd48f76",
   "metadata": {},
   "outputs": [],
   "source": [
    "# split data\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b778bd16",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import GradientBoostingRegressor\n",
    "from sklearn.metrics import r2_score\t\n",
    "\n",
    "\n",
    "# Fit the model on the filtered Train set and print summary\n",
    "model_4 = GradientBoostingRegressor(n_estimators=1_000, random_state=42)\n",
    "fitted_4 = model_4.fit(X_train, y_train)\n",
    "y_pred_gb_0 = model_4.predict(X_test)\n",
    "\n",
    "# Calculate MAE for the test set\n",
    "mae_gb = mean_absolute_error(y_true=y_test, y_pred=y_pred_gb_0)\n",
    "print(f\"GB MAE: ${mae_gb:,.2f}\")\n",
    "\n",
    "r2_gb = r2_score(y_true=y_test, y_pred=y_pred_gb_0)\n",
    "print(f\"GB R-squared: {r2_gb:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dad1864c",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "We can see that the model performs better than the OLS model. But we still need to check for heteroskedasticity.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f340975a",
   "metadata": {},
   "outputs": [],
   "source": [
    "resid_gb_0 = y_train - fitted_4.predict(X_train)\n",
    "\n",
    "# plot residuals vs. predicted values\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "sns.scatterplot(x=y_train, y=resid_gb_0, color='royalblue', edgecolor='black', alpha=0.8)\n",
    "plt.axhline(y=0, color='r', linestyle='--', linewidth=1)\n",
    "plt.title(\"Residuals vs. Predicted Values\\n (On Test Set)\")\n",
    "plt.xlabel(\"Predicted Values\")\n",
    "plt.xticks(rotation=0)\n",
    "plt.ylabel(\"Residuals\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ee84d2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform Breusch-Pagan test\n",
    "bp_test = het_breuschpagan(resid_gb_0, sm.add_constant(X_train)) # type: ignore\n",
    "print(\"Breusch-Pagan test p-value:\", bp_test[1])\n",
    "\n",
    "# Check if p-value is less than the significance level\n",
    "alpha = 0.05\n",
    "if bp_test[1] < alpha:\n",
    "    print(f\"Reject the null hypothesis: Heteroskedasticity is present.\")\n",
    "else:\n",
    "    print(f\"Fail to reject the null hypothesis: Heteroskedasticity is not present.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "073b0f6a",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "Unfortunately, we still observe heteroskedasticity in the residuals both using the `Breusch-Pagan` test and the `Visual Inspection` technique.\n",
    "\n",
    "### **Fitting the model on Filtered Data**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8624a91",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit the model on the filtered Train set and print summary\n",
    "model_5 = GradientBoostingRegressor(n_estimators=1_000, random_state=42)\n",
    "\n",
    "fitted_5 = model_5.fit(filtered_X_train.drop('const', axis=1), filtered_y_train)  # type: ignore\n",
    "y_pred_gb_1 = model_5.predict(X_test)\n",
    "\n",
    "# Calculate MAE for the test set\n",
    "mae_gb_2 = mean_absolute_error(y_true=y_test, y_pred=y_pred_gb_1)\n",
    "print(f\"GB MAE: ${mae_gb_2:,.2f}\")\n",
    "\n",
    "r2_gb_2= r2_score(y_true=y_test, y_pred=y_pred_gb_1)\n",
    "print(f\"GB R-squared: {r2_gb_2:.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73dcc537",
   "metadata": {},
   "outputs": [],
   "source": [
    "resid_gb_1 = y_train - fitted_5.predict(X_train)\n",
    "\n",
    "# Perform Breusch-Pagan test\n",
    "bp_test = het_breuschpagan(resid_gb_1, sm.add_constant(X_train)) # type: ignore\n",
    "print(\"Breusch-Pagan test p-value:\", bp_test[1])\n",
    "\n",
    "# Check if p-value is less than the significance level\n",
    "alpha = 0.05\n",
    "if bp_test[1] < alpha:\n",
    "    print(f\"Reject the null hypothesis: Heteroskedasticity is present.\")\n",
    "else:\n",
    "    print(f\"Fail to reject the null hypothesis: Heteroskedasticity is not present.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9d69fb2",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "The model still has heteroskedasticity.\n",
    "### **Fitting the model on Transformed & Scaled Data**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21dde887",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_sc.drop('const', axis=1, inplace=True)  # type: ignore\n",
    "X_test_sc.drop('const', axis=1, inplace=True)  # type: ignore"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3b56cce",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit the model on the filtered Train set and print summary\n",
    "model_6 = GradientBoostingRegressor(n_estimators=1_000, random_state=42)\n",
    "\n",
    "fitted_6 = model_6.fit(X_train_sc, y_train_sc)  # type: ignore\n",
    "y_pred_gb_2 = model_6.predict(X_test_sc)\n",
    "y_pred__gb_2_inv = inv_boxcox(y_pred_gb_2, best_lambda)\n",
    "\n",
    "# Calculate MAE for the test set\n",
    "mae_gb_3 = mean_absolute_error(y_true=y_test, y_pred=y_pred__gb_2_inv)\n",
    "print(f\"GB MAE: ${mae_gb_3:,.2f}\")\n",
    "\n",
    "r2_gb_3= r2_score(y_true=y_test, y_pred=y_pred__gb_2_inv)\n",
    "print(f\"GB R-squared: {r2_gb_3:.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc636c50",
   "metadata": {},
   "outputs": [],
   "source": [
    "resid_gb_2 = y_train_sc - fitted_6.predict(X_train_sc)\n",
    "\n",
    "# Perform Breusch-Pagan test\n",
    "bp_test = het_breuschpagan(resid_gb_2, sm.add_constant(X_train_sc)) # type: ignore\n",
    "print(\"Breusch-Pagan test p-value:\", bp_test[1])\n",
    "\n",
    "# Check if p-value is less than the significance level\n",
    "alpha = 0.05\n",
    "if bp_test[1] < alpha:\n",
    "    print(f\"Reject the null hypothesis: Heteroskedasticity is present.\")\n",
    "else:\n",
    "    print(f\"Fail to reject the null hypothesis: Heteroskedasticity is not present.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "884162b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot residuals vs. predicted values\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "sns.scatterplot(x=y_train_sc, y=resid_gb_2, color='royalblue', edgecolor='black', alpha=0.8)  # type: ignore\n",
    "plt.axhline(y=0, color='r', linestyle='--', linewidth=1)\n",
    "plt.title(\"Residuals vs. Predicted Values\\n (On Test Set)\")\n",
    "plt.xlabel(\"Predicted Values\")\n",
    "plt.xticks(rotation=0)\n",
    "plt.ylabel(\"Residuals\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "206df4a3",
   "metadata": {},
   "source": [
    "<div style=\"font-family: 'Times New Roman', serif; w\">\n",
    "\n",
    "We can observe that the model has improved and heteroskedasticity is no longer present. \n",
    "\n",
    "Can we do better?\n",
    "\n",
    "#### **Fine Tuning the Model: Leverage & Residuals Threshold**\n",
    "\n",
    "We wil build a function that takes the influence_df as input and identify high leverage and high residual points. These points will be removed from the Train set with the constraint of high `Breusch-Pagan test p-value` and lower `MAE`.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fca128d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_leverage_and_residuals_threshold(influence_df, threshold_leverage, threshold_residual):\n",
    "    df_results = pd.DataFrame(columns=['high_leverage_threshold', 'high_residual_threshold', 'Breusch-Pagan test p-value', 'MAE', 'R2'])\t\n",
    "    \n",
    "    for i in threshold_leverage:\n",
    "        # append threshold_leverage to df_results\n",
    "        for j in threshold_residual:  \n",
    "\n",
    "            # append threshold_residual to df_results\n",
    "            high_leverage_threshold = i\n",
    "            high_residual_threshold = j\n",
    "\n",
    "            # Identify high leverage and high residual points\n",
    "            high_leverage_points = influence_df.cooks_d > high_leverage_threshold\n",
    "            high_residual_points = np.abs(influence_df.standard_resid) > high_residual_threshold\n",
    "\n",
    "            # Combine both conditions to find points with either high leverage or high residuals\n",
    "            influential_points = high_leverage_points & high_residual_points \n",
    "\n",
    "            # Remove influential points from the Train set\n",
    "            filtered_X_train = X_train[~influential_points]\n",
    "            filtered_y_train = y_train[~influential_points]\n",
    "\n",
    "            # Apply Box-Cox transformation to `price_usd` & `covered_surface_in_m2`\n",
    "            X_train_sc = filtered_X_train.copy()\n",
    "            y_train_sc = filtered_y_train.copy()\n",
    "\n",
    "            X_test_sc = X_test.copy()\n",
    "            y_test_sc = y_test.copy()\n",
    "\n",
    "            best_lambda = test_lambda_df.lambda_value[0]\n",
    "            X_train_sc.surface_covered_in_m2 = sps.boxcox(X_train_sc.surface_covered_in_m2, best_lambda)\n",
    "            y_train_sc = sps.boxcox(y_train_sc, best_lambda)\n",
    "\n",
    "            X_test_sc.surface_covered_in_m2 = sps.boxcox(X_test_sc.surface_covered_in_m2, best_lambda)\n",
    "            y_test_sc = sps.boxcox(y_test_sc, best_lambda)\n",
    "\n",
    "            # Scale Data\n",
    "            std_scaler = StandardScaler(with_mean=True, with_std=True, copy=True)\n",
    "\n",
    "            X_train_sc = pd.DataFrame(std_scaler.fit_transform(X_train_sc), columns=X_train.columns, index=filtered_X_train.index)\n",
    "            X_test_sc = pd.DataFrame(std_scaler.transform(X_test_sc), columns=X_test.columns, index=X_test.index)  # type: ignore\n",
    "\n",
    "            # Fit the model on the filtered Train set and print summary\n",
    "            model = GradientBoostingRegressor(n_estimators=1_000, random_state=42)\n",
    "\n",
    "            fitted= model.fit(X_train_sc, y_train_sc) # type: ignore\n",
    "            y_pred_gb= model.predict(X_test_sc)\n",
    "            y_pred__gb_inv = inv_boxcox(y_pred_gb, best_lambda)\n",
    "\n",
    "            # Calculate MAE for the test set\n",
    "            mae_gb = mean_absolute_error(y_true=y_test, y_pred=y_pred__gb_inv)\n",
    "            r2_gb= r2_score(y_true=y_test, y_pred=y_pred__gb_inv)\n",
    "\n",
    "            resid_train = y_train_sc - fitted.predict(X_train_sc) \n",
    "\n",
    "            bp_test = het_breuschpagan(resid_train, sm.add_constant(X_train_sc))[1] # type: ignore\n",
    "\n",
    "            new_row = pd.DataFrame({\n",
    "                'high_leverage_threshold': i,\n",
    "                'high_residual_threshold': j,\n",
    "                'Breusch-Pagan test p-value': bp_test,\n",
    "                'MAE': mae_gb,\n",
    "                'R2': r2_gb}, index=[0]) # type: ignore\n",
    "            df_results = pd.concat([df_results, new_row], ignore_index=True)\n",
    "\n",
    "    return df_results.sort_values(by='MAE', ascending=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eda9f7e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold_leverage = np.arange(0.0003, 0.0004, 0.00005)\n",
    "threshold_residual = np.arange(3, 3.6, 0.5)\n",
    "\n",
    "finale_results = test_leverage_and_residuals_threshold(influence_df, threshold_leverage, threshold_residual)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6494665d",
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.options.display.float_format = '{:.4f}'.format\n",
    "finale_results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a02d3bac",
   "metadata": {},
   "source": [
    "<div style=\"font-family: Times New Roman, serif;\">\n",
    "\n",
    "Our final result are: the best `high_leverage_threshold` is **0.0003** and the best `high_residual_threshold` is **3**. We will proceed with these two values as they will give us the best results, resulting in a MAE of **$23,809.9879** and an R2 of **0.7863**.\n",
    "\n",
    "The best part is that we successfully handled the heteroskedasticity issue.\n",
    "\n",
    "Next, we will store our pretrained model in a pickle file.\n",
    "\n",
    "## **Store Model**\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fe0484f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Identify high leverage and high residual points\n",
    "high_leverage_points = influence_df.cooks_d > 0.0003\t\n",
    "high_residual_points = np.abs(influence_df.standard_resid) > 3\n",
    "\n",
    "# Combine both conditions to find points with either high leverage or high residuals\n",
    "influential_points = high_leverage_points & high_residual_points \n",
    "\n",
    "# Remove influential points from the Train set\n",
    "filtered_X_train = X_train[~influential_points]\n",
    "filtered_y_train = y_train[~influential_points]\n",
    "\n",
    "# Apply Box-Cox transformation to `price_usd` & `covered_surface_in_m2`\n",
    "X_train_sc = filtered_X_train.copy()\n",
    "y_train_sc = filtered_y_train.copy()\n",
    "\n",
    "X_train_sc.surface_covered_in_m2 = sps.boxcox(X_train_sc.surface_covered_in_m2, best_lambda)\n",
    "y_train_sc = pd.Series(sps.boxcox(y_train_sc, best_lambda), index=X_train_sc.index, name='price_aprox_usd')\n",
    "\n",
    "# Scale Data\n",
    "std_scaler = StandardScaler(with_mean=True, with_std=True, copy=True)\n",
    "\n",
    "X_train_sc = pd.DataFrame(std_scaler.fit_transform(X_train_sc), columns=X_train.columns, index=filtered_X_train.index)\n",
    "\n",
    "# Fit the model on the filtered Train set and print summary\n",
    "model = GradientBoostingRegressor(n_estimators=1_000, random_state=42)\n",
    "\n",
    "fitted= model.fit(X_train_sc, y_train_sc) # type: ignore"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6486e252",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Store our model in a pickle file\n",
    "import pickle\n",
    "\n",
    "pickle.dump(fitted, open('final_model.pkl', 'wb'))\n",
    "pickle.dump(std_scaler, open('standard_scaler.pkl', 'wb'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17c4cd14",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load our model from the pickle file\n",
    "model = pickle.load(open('final_model.pkl', 'rb'))\n",
    "std_scaler = pickle.load(open('standard_scaler.pkl', 'rb'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82da0337",
   "metadata": {},
   "source": [
    "<div style=\"font-family: Times New Roman, serif;\">\n",
    "\n",
    "## **Deployment**\n",
    "\n",
    "Our project is a Dash application that provides a real estate price prediction tool. It uses the Dash framework to create a web application with a user interface. The layout of the app is defined using HTML and CSS, and it includes various components such as input fields, buttons, and a map.\n",
    "\n",
    "The app allows users to input the covered surface, latitude, and longitude of a property. When the user clicks the \"Predict Price\" button, the app retrieves the address based on the latitude and longitude using a function called `get_address()`. It then checks if the city is in Argentina and if it is in the Autonomous City of Buenos Aires. If not, it displays an error message.\n",
    "\n",
    "If the city is valid, the app prepares a feature vector using the input values and performs a price prediction using a pre-trained model. The predicted price is then displayed on the app.\n",
    "\n",
    "The app also includes a map component that allows users to click on a location to get the latitude and longitude coordinates. The coordinates are displayed on the app and can be used as input for the price prediction.\n",
    "\n",
    "Additionally, the app provides information about the municipality where the property is located. It retrieves data from our data and displays statistics about the price distribution in that municipality.\n",
    "\n",
    "The `get_address()` function retrieves the address based on the latitude and longitude using the `geopy` library.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c753cc7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from geopy.geocoders import Nominatim  # type: ignore\n",
    "\n",
    "def get_address(latitude, longitude):\n",
    "    geolocator = Nominatim(user_agent=\"coordinateconverter\")\n",
    "    try:\n",
    "        location = geolocator.reverse((latitude, longitude), exactly_one=True, timeout=10, addressdetails=True, language='en', zoom=16)  # type: ignore\n",
    "        return location.address  # type: ignore\n",
    "    except:\n",
    "        return \"Not Found!\"\n",
    "\n",
    "# Example usage\n",
    "latitude = -34.617768\t\t\n",
    "longitude = -58.566394\t\n",
    "address = get_address(latitude, longitude)\n",
    "print(address)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2326ee08",
   "metadata": {},
   "outputs": [],
   "source": [
    "from dash import Dash, html, dcc, Input, Output, State, dash_table\n",
    "from dash_extensions.javascript import assign\n",
    "import dash_leaflet as dl\n",
    "\n",
    "# Initialize the Dash app\n",
    "app = Dash(__name__)\n",
    "\n",
    "# JavaScript code to capture click event and set coordinates as property\n",
    "eventHandlers = dict(\n",
    "    click=assign(\"function(e, ctx){ctx.setProps({data: {lat: e.latlng.lat, lng: e.latlng.lng}})}\")\n",
    ")\n",
    "\n",
    "# Define the layout of the app\n",
    "app.layout = html.Div([\n",
    "    html.Div([\n",
    "        html.H1(\"Real Estate Price Prediction\", style={'text-align': 'center'}),\n",
    "    ], style={'font-weight': 'bold', 'font-family':'Calibri','color':'#FD8916', 'background-color': '#17616E', 'width': '60%', 'padding':'1px','text-align': 'center', 'margin': '0 auto', 'border': '0px solid white', 'border-radius': '10px', 'margin-bottom': '10px'}),\n",
    "\n",
    "    html.Div([\n",
    "        html.Div([\n",
    "            html.Div([\n",
    "                html.Div([\n",
    "                    dcc.Input(id='covered_surface', type='number', placeholder='Surface Covered in m\\u00b2',\n",
    "                            style={'margin': '10px', 'padding': '10px', 'border-radius': '8px', 'border': '1px solid #FD8916'}),\n",
    "\n",
    "                    dcc.Input(id='latitude', type='number', placeholder='Latitude',\n",
    "                            style={'margin': '10px', 'padding': '10px', 'border-radius': '8px', 'border': '1px solid #FD8916'}),\n",
    "\n",
    "                    dcc.Input(id='longitude', type='number', placeholder='Longitude',\n",
    "                            style={'margin': '10px', 'padding': '10px', 'border-radius': '8px', 'border': '1px solid #FD8916'}),\n",
    "\n",
    "                ], style={'width': '100%', 'border': '0px solid #007BFF', 'border-radius': '10px'}),\n",
    "\n",
    "                # Define the button\n",
    "                html.Div([\n",
    "                    html.Button('Predict Price', id='predict_button', n_clicks=0,\n",
    "                            style={'font-weight': 'bold', 'margin': '10px', 'padding': '10px', 'background-color': '#FD8916', 'color': '#173b61', 'border-radius': '8px', 'border': '1px solid #173b61', 'float': 'right'}),\n",
    "                ], style={'width': '100%', 'border': '0px solid #007BFF', 'border-radius': '10px', 'display': 'flex', 'justify-content': 'right'}),\n",
    "            ], style={'border': '2px solid #17616E', 'border-radius': '10px'}),\n",
    "\n",
    "            html.Div(id='prediction_output', style={'width': '100%', 'height': 'auto', 'margin-top': '10px','border': '2px solid #17616E', 'border-radius': '10px', 'padding': '0px', 'display':'none'},\n",
    "                    children=[\n",
    "                        html.Div([\n",
    "                            html.Div(\"Predicted Price:\", style={'display': 'inline-block', 'font-weight': 'bold', 'font-family':'Calibri'}),\n",
    "                            html.Div(id='predicted_price', style={'font-weight': 'bold', 'font-family':'Calibri', 'display': 'inline-block', 'margin-left': '20px', 'color':'#FD8916'}),\n",
    "                        ], style={'padding': '10px', 'margin-bottom': '10px'}),\n",
    "                        html.Div([\n",
    "                            html.Div(\"Apartment Address.:\", style={'display': 'inline-block', 'font-weight': 'bold', 'font-family':'Calibri'}),\n",
    "                            html.Div(id='address', style={'font-weight': 'bold', 'font-family':'Calibri', 'color': '#17616E', 'display': 'inline-block', 'margin-left': '10px'}),\n",
    "                        ], style={'padding': '10px', 'margin-bottom': '10px'}),\n",
    "                        html.Div([\n",
    "                            html.Div(\"Municipality Statistics:\", style={'display': 'inline-block', 'font-weight': 'bold', 'font-family':'Calibri'}),\n",
    "                            html.Div(id='municipality_name', style={'font-weight': 'bold', 'font-family':'Calibri','color':'#FD8916', 'margin': '5px', 'display': 'flex', 'justify-content': 'center'}),\n",
    "                            html.Div(id='municipality_info', style={'font-weight': 'bold', 'font-family':'Calibri', 'color': '#17616E', 'margin': '10px', 'display': 'flex', 'justify-content': 'center'}),\n",
    "                        ], style={'padding': '10px', 'margin-bottom': '10px'}),\n",
    "                    ]\n",
    "            ),\n",
    "        ], style={'width': '40%', 'padding': '20px', 'border': '2px solid #17616E', 'border-radius': '10px', 'margin-right': '10px', 'background-color': 'white'}),\n",
    "        \n",
    "        html.Div([\n",
    "            dl.Map(\n",
    "                children=[\n",
    "                    dl.TileLayer(),\n",
    "                    dl.Marker(id='marker', position=[33.6086, -7.6327])],\n",
    "                eventHandlers=eventHandlers,\n",
    "                style={'height': '60vh', 'width': '100%'},  # Set height to 80% of the viewport height\n",
    "                center=[-34.6037, -58.3816],  # Center around Buenos Aires\n",
    "                zoom=12,  # Set zoom level to view the city\n",
    "                id='map'\n",
    "            ),\n",
    "            html.Div(id='coordinates', style={'marginTop': '20px', 'font-weight': 'bold', 'font-family':'Calibri', 'fontSize': '20px', 'color': 'black', 'backgroundColor': 'white', 'border-radius': '10px', 'padding': '10px', 'text-align': 'center'})\n",
    "        ], style={'width': '60%', 'padding': '20px', 'border': '0px solid #FAF3DD', 'border-radius': '10px', 'background-color': '#173b61'}\n",
    "        ),\n",
    "    ], style={'display': 'flex', 'justify-content': 'space-between'})\n",
    "], style={'background-color': 'white', 'max-width': '1200px', 'margin': '0 auto', 'backgroundColor': '#FFEBD0', 'padding': '20px', 'border': '2px solid #17616E', 'border-radius': '10px'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc7ce361",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the callback to show the prediction output when the button is clicked\n",
    "@app.callback(\n",
    "    Output('prediction_output', 'style'),\n",
    "    Input('predict_button', 'n_clicks'),\n",
    "    prevent_initial_call=False \n",
    ")\n",
    "def show_prediction_output(n_clicks):\n",
    "    if n_clicks > 0:\n",
    "        return {'display': 'block', 'width': '100%', 'height': 'auto', 'margin-top': '10px', 'border': '2px solid #17616E', 'border-radius': '10px', 'padding': '0px'}\n",
    "    return {'display': 'none'}\n",
    "\n",
    "# Callback to update coordinates based on map click\n",
    "@app.callback(\n",
    "    [Output('latitude', 'value'),\n",
    "     Output('longitude', 'value'),\n",
    "     Output('coordinates', 'children'),\n",
    "     Output('marker', 'position')],\n",
    "    Input('map', 'data')\n",
    ")\n",
    "def update_coordinates(data):\n",
    "    if data:\n",
    "        return data['lat'], data['lng'], f\"Coordinates: Latitude {data['lat']:.4f}, Longitude {data['lng']:.4f}\", [data['lat'], data['lng']]\n",
    "    return None, None, \"Click on the map to get coordinates\", [33.6086, -7.6327]\n",
    "\n",
    "# Separate callback to update address\n",
    "@app.callback(\n",
    "    Output('address', 'children'),\n",
    "    [Input('predict_button', 'n_clicks')],\n",
    "    [State('latitude', 'value'),\n",
    "     State('longitude', 'value')]\n",
    ")\n",
    "def update_address(n_clicks, latitude, longitude):\n",
    "    if n_clicks > 0 and all(v is not None for v in [latitude, longitude]):\n",
    "        # Get address\n",
    "        address = get_address(latitude, longitude)\n",
    "        return address\n",
    "    return ''\n",
    "\n",
    "# Separate callback to update predicted price\n",
    "@app.callback(\n",
    "    Output('predicted_price', 'children'),\n",
    "    [Input('predict_button', 'n_clicks')],\n",
    "    [State('covered_surface', 'value'),\n",
    "     State('latitude', 'value'),\n",
    "     State('longitude', 'value')]\n",
    ")\n",
    "def update_predicted_price(n_clicks, covered_surface, latitude, longitude):\n",
    "    if n_clicks > 0 and all(v is not None for v in [covered_surface, latitude, longitude]):\n",
    "\n",
    "        # check if city is in Argentina by using the get_address function\n",
    "        try:\n",
    "            city = [element.strip() for element in get_address(latitude, longitude).split(\",\")][-3]\n",
    "        except:\n",
    "            return 'So you want to buy an apartment in the sea... Very funny 😄'\n",
    "\n",
    "        if city != 'Autonomous City of Buenos Aires':\n",
    "            return '❌ The city is not in the Autonomous City of Buenos Aires! Please select another location.'\n",
    "\n",
    "        # Prepare the feature vector\n",
    "        features = pd.DataFrame([[latitude, longitude, covered_surface]],  # the 'price_aprox_usd' column is set to 0, it won't affect the prediction, it's just a placeholder\n",
    "                                columns=['latitude', 'longitude', 'surface_covered_in_m2'])\n",
    "\n",
    "        # Transform the feature `surface_covered_in_m2` using the 'box-cox' transformation\n",
    "        features.surface_covered_in_m2 = sps.boxcox(features.surface_covered_in_m2, best_lambda)\n",
    "\n",
    "        # Scale the feature vector\n",
    "        features_scaled = std_scaler.transform(features)      \n",
    "        features_df = pd.DataFrame(features_scaled, columns=features.columns)  # type: ignore\n",
    "\n",
    "        # Predict Apartment Price \n",
    "        prediction = model.predict(features_df)\n",
    "\n",
    "        # Inverse Box-Cox transformation\n",
    "        prediction_orginal_scale = inv_boxcox(prediction, best_lambda)[0]\n",
    "\n",
    "        # Return the predicted price\n",
    "        return f'${prediction_orginal_scale:,.0f}'\n",
    "    return ''\n",
    "\n",
    "# Separate callback to update municipality information\n",
    "@app.callback(\n",
    "    [Output('municipality_name', 'children'),\n",
    "    Output('municipality_info', 'children')],\n",
    "    [Input('predict_button', 'n_clicks')],\n",
    "    [State('latitude', 'value'),\n",
    "     State('longitude', 'value')]\n",
    ")\n",
    "def update_municipality_info(n_clicks, latitude, longitude):\n",
    "    if n_clicks > 0 and all(v is not None for v in [latitude, longitude]):\n",
    "        try:\n",
    "            # Get address\n",
    "            address = get_address(latitude, longitude)\n",
    "\n",
    "            # Filter data by municipality\n",
    "            municipality = [element.strip().lower() for element in address.split(\",\")]\n",
    "            info_municipality = df_full[(df_full.municipality.str.lower().str.contains(municipality[1]))\n",
    "                                        | (df_full.municipality.str.lower().str.contains(municipality[2]))\n",
    "                                        | (df_full.municipality.str.lower().str.contains(municipality[3]))\n",
    "            ]\n",
    "        except Exception:\n",
    "            return '', f\"🔱 Poseidon has put up a ‘No Humans Allowed’ sign for the sea apartments. Apparently, he's not into ocean-side neighbors who don't appreciate underwater real estate!\"\n",
    "        \n",
    "        if info_municipality.shape[0] == 0:\n",
    "            return '', f\"🤷🏻‍♂️ It seems that we do not have any information about this location\"\n",
    "        else:\n",
    "            # Extract municipality name\n",
    "            municipality_name = info_municipality.municipality.unique()[0]\n",
    "\n",
    "            # Create table\n",
    "            municipality_df = pd.DataFrame(info_municipality.price_aprox_usd.describe())[1:]  # Exclude 'count'\n",
    "\n",
    "            # Create 'Measure' column\n",
    "            municipality_df['Measure'] = 0\n",
    "\n",
    "            # Order columns & rename `price_usd` to `Value`\n",
    "            municipality_df = municipality_df[['Measure', 'price_aprox_usd']].rename(columns={'price_aprox_usd': 'Value'})  # type: ignore \n",
    "\n",
    "            # Format table\n",
    "            municipality_df.Value = municipality_df.Value.apply(lambda x: f\"${x:,.0f}\" if isinstance(x, (int, float)) else x)  # type: ignore\n",
    "\n",
    "            # Rename rows to proper headings\n",
    "            municipality_df['Measure'] = ['Average Price (USD)', 'Price Standard Deviation (USD)', 'Minimum Price (USD)', '25th Percentile Price (USD)', \n",
    "                                    'Median Price (USD)', '75th Percentile Price (USD)','Maximum Price (USD)']\n",
    "\n",
    "            # Return the table        \n",
    "            return municipality_name, dash_table.DataTable(\n",
    "                municipality_df.to_dict('records'),[{'name': 'Measure', 'id': 'Measure'}, {'name': 'Value', 'id': 'Value'}],\n",
    "                style_cell={'textAlign': 'left', 'font-family': 'Calibri', 'padding': '10px', 'border': '1px solid #17616E'},\n",
    "                style_header={'backgroundColor': '#17616E', 'color': 'white', 'fontWeight': 'bold', 'textAlign': 'center'},\n",
    "                style_table={'width': '100%', 'margin': '0 auto'}\n",
    "                )\n",
    "        \n",
    "    return '', html.Div('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7fc5cb8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the app\n",
    "if __name__ == '__main__':\n",
    "    app.run_server(debug=True, port=8080)\n",
    "\n",
    "# Open http://127.0.0.1:8080/"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "405fb582",
   "metadata": {},
   "source": [
    "<div style=\"width: 90%; font-family: Times New Roman, serif; text-justify: inter-word; margin: 0 auto; font-size: 20px;\">\n",
    "\n",
    "## **Conclusion**\n",
    "\n",
    "In conclusion, our project **successfully** addressed key challenges in predicting real estate prices by effectively handling **heteroskedasticity** and optimizing model performance. Through rigorous analysis, we determined that the optimal `high_leverage_threshold` is **0.0003** and the best `high_residual_threshold` is **3**. Adopting these thresholds led to significant improvements in our model, achieving a **Mean Absolute Error (MAE)** of ***$23,809.9879*** and an **R²** of ***0.7863***.\n",
    "\n",
    "Notably, we surpassed the baseline MAE by **60%**, demonstrating the effectiveness of our approach. Extensive **feature engineering** and the application of advanced techniques enabled us to identify and address outliers and influential points, ultimately refining our model.\n",
    "\n",
    "Our iterative process involved experimenting with various versions of the **Ordinary Least Squares (OLS)** model, each iteration contributing to the overall enhancement of model performance. The final choice of a **gradient boosting** model was informed by its superior performance and ability to address heteroskedasticity, ensuring the most accurate predictions and robust results.\n",
    "\n",
    "Overall, our approach and methodology provided a well-rounded solution to the complexities of real estate price prediction, showcasing our **capability** to leverage data science techniques to drive impactful insights.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92250fdf",
   "metadata": {},
   "source": [
    "<div style=\"width: 90%; font-family: Times New Roman, serif; text-justify: inter-word; margin: 0 auto; font-size: 20px;\">\n",
    "\n",
    "## **References**\n",
    "\n",
    "\n",
    "1. **Data Sources**:\n",
    "   - [Buenos Aires Real Estate Data]. (2016). Retrieved from [Properati.com](https://properati.com/)\n",
    "\n",
    "2. **Technical Papers and Articles**:\n",
    "   - Christian Thieme. (2021). Identifying Outliers in Linear Regression — Cook’s Distance. [Medium.com](https://towardsdatascience.com/identifying-outliers-in-linear-regression-cooks-distance-9e212e9136a).\n",
    "   - Eric Hansen. (2021). Linear Regression and Assumption Validity.[Medium.com](https://ericthansen.medium.com/linear-regression-and-assumption-validity-94712f714ea4).\n",
    "   - Hitesh Taneja. (2023). In-Depth Guide to Cook’s Distance. [Medium.com](https://medium.com/@hiteshtaneja/in-depth-guide-to-cooks-distance-40697cda2c9d).\n",
    "\n",
    "3. **Books and Tutorials**:\n",
    "   - Aman Kharwal. (2023). Build An End-to-End Machine Learning Model. [thecleverprogrammer.com](https://thecleverprogrammer.com/2023/1218build-an-end-to-end-machine-learning-model/).\n",
    "\n",
    "4. **Official Documentations**:\n",
    "   - Events. [dash-leaflet.com](https://www.dash-leaflet.com/docs/events)\n",
    "   - Dash DataTable. [Dash.Plotly.com](https://dash.plotly.com/datatable)\n",
    "   - Performance. [Dash.Plotly.com](https://dash.plotly.com/performance)\n",
    "\n",
    "\n",
    "5. **Libraries and Tools**:\n",
    "   - Scikit-learn.\n",
    "   - Statsmodels.\n",
    "   - Pandas.\n",
    "   - NumPy.\n",
    "   - Matplotlib.\n",
    "   - Seaborn.\n",
    "   - Plotly.\n",
    "   - Dash.\n",
    "\n",
    "## **About the Author**\n",
    "\n",
    "**Youssef Laouina** is a certified data scientist with expertise in Python, machine learning, and data analysis. He has a background in economics and has worked on projects addressing real-world problems through data-driven solutions.\n",
    "\n",
    "Connect with Youssef:\n",
    "- **Email**: [laouina.yusuf@gmail.com](mailto:laouina.yusuf@gmail.com)\n",
    "- **LinkedIn**: [Youssef Laouina](https://www.linkedin.com/in/youssef-laouina)\n",
    "- **GitHub**: [Youssef Laouina](https://github.com/youssef-laouina)\n",
    "\n",
    "   </div>\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.2"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}