causal_inference_with_vertex_ai_automl_forecasting.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Causal Inference with Vertex AI AutoML Forecasting.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "oDb0nyJ5j4Q0"
      },
      "source": [
        "# Copyright 2021 Google LLC\n",
        "#\n",
        "# Licensed under the Apache License, Version 2.0 (the \"License\");\n",
        "# you may not use this file except in compliance with the License.\n",
        "# You may obtain a copy of the License at\n",
        "#\n",
        "#     https://www.apache.org/licenses/LICENSE-2.0\n",
        "#\n",
        "# Unless required by applicable law or agreed to in writing, software\n",
        "# distributed under the License is distributed on an \"AS IS\" BASIS,\n",
        "# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
        "# See the License for the specific language governing permissions and\n",
        "# limitations under the License."
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1W3qZ2qHj7Qs"
      },
      "source": [
        "<table align=\"left\">\n",
        "  <td>\n",
        "    <a href=\"https://console.cloud.google.com/ai-platform/notebooks/deploy-notebook?name=Vertex%20AI%20Pipeline&download_url=https%3A%2F%2Fraw.githubusercontent.com%2Fkweinmeister%2Fnotebooks%2Fmaster%2Fcausal_inference_with_vertex_ai_automl_forecasting.ipynb\">\n",
        "      <img src=\"https://cloud.google.com/images/products/ai/ai-solutions-icon.svg\" alt=\"Google Cloud Notebooks\"> Open in GCP Notebooks\n",
        "    </a>\n",
        "  </td> \n",
        "  <td>\n",
        "    <a href=\"https://colab.research.google.com/github/kweinmeister/notebooks/blob/master/causal_inference_with_vertex_ai_automl_forecasting.ipynb\">\n",
        "      <img src=\"https://cloud.google.com/ml-engine/images/colab-logo-32px.png\" alt=\"Colab logo\"> Open in Colab\n",
        "    </a>\n",
        "  </td>\n",
        "  <td>\n",
        "    <a href=\"https://github.com/kweinmeister/notebooks/blob/master/causal_inference_with_vertex_ai_automl_forecasting.ipynb\">\n",
        "        <img src=\"https://cloud.google.com/ml-engine/images/github-logo-32px.png\" alt=\"GitHub logo\">\n",
        "      View on GitHub\n",
        "    </a>\n",
        "  </td>\n",
        "</table>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kmvZlc-i6rG6"
      },
      "source": [
        "# Causal Inference with Vertex AI AutoML Forecasting"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tgV4bZmG654u"
      },
      "source": [
        "The goal of this notebook is to demonstrate how to estimate the effect of a treatment. A treatment, or intervention, is an action intended to change an outcome. For example, you could estimate the effect of a marketing campaign or promotion on sales.\n",
        "\n",
        "We'll look at answering the question, \"How did Brexit impact exchange rates between the British Pound and US Dollar?\" We'll use time-series data from [FRED](https://fred.stlouisfed.org).\n",
        "\n",
        "This notebook will use the tabular data forecasting service from [Vertex AI](https://cloud.google.com/vertex-ai):\n",
        "* We will train a time-series forecasting model on the British Pound and a related time-series (Euro) to generate a counterfactual.\n",
        "* The counterfactual time series aims to answer the question \"What would the value of the Pound be, had this event not happened?\"\n",
        "* We will then predict the exchange rate during the 30 day period after the Brexit vote on June 23, 2016.\n",
        "* The difference between the actual and counterfactual time series gives us an estimate of the vote.\n",
        "\n",
        "Finally, we'll use [tfcausalimpact](https://github.com/WillianFuks/tfcausalimpact) to compare our results to a statistical approach. tfcausalimpact is a Python port of the R-based [CausalImpact](https://google.github.io/CausalImpact/), using [TensorFlow Probability](https://www.tensorflow.org/probability).\n",
        "\n",
        "Citation:\n",
        "* Board of Governors of the Federal Reserve System (US), retrieved from FRED, Federal Reserve Bank of St. Louis:\n",
        "  * [UK / US Exchange Rate - DEXUSUK](https://fred.stlouisfed.org/series/DEXUSUK)\n",
        "  * [Euro / US Exchange Rate - DEXUSEU](https://fred.stlouisfed.org/series/DEXUSEU)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dwrxP2fPszac"
      },
      "source": [
        "## Imports"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "f-Jwo7A0wSY0"
      },
      "source": [
        "# Install google-cloud-aiplatform (and other packages to match)\n",
        "\n",
        "!pip3 install -U --quiet google-cloud-aiplatform google-cloud-bigquery grpcio pandas-gbq gcsfs tfcausalimpact folium==0.2.1 imgaug==0.2.6 tensorflow==2.5\n",
        "\n",
        "# Restart runtime\n",
        "\n",
        "import IPython\n",
        "app = IPython.Application.instance() \n",
        "app.kernel.do_shutdown(True)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BDkQ8z-0N6ht"
      },
      "source": [
        "# Import packages\n",
        "\n",
        "from causalimpact import CausalImpact\n",
        "from datetime import timedelta\n",
        "from scipy import stats\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import sys"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PGZJA4las6V1"
      },
      "source": [
        "## Constants"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "x9C2vKqYu3Is"
      },
      "source": [
        "# Set GCP project constants\n",
        "\n",
        "PROJECT_ID = 'YOUR-PROJECT-ID' # Your project ID\n",
        "REGION = 'us-central1' # Your project region\n",
        "STORAGE_PATH = 'gs://YOUR-REGIONAL-BUCKET/data/exchange_rates' # A bucket and path to store various CSV files used (training, forecast). Bucket must already exist."
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ChRDzsmiWc0Q"
      },
      "source": [
        "# Set forecasting constants (no need to change these)\n",
        "\n",
        "# Training data prior to event\n",
        "PRE_PERIOD=[pd.to_datetime('2012-01-03'), pd.to_datetime('2016-06-23')]\n",
        "\n",
        "# Duration of period after the event to analyze\n",
        "POST_PERIOD_LENGTH=28\n",
        "POST_PERIOD=[pd.to_datetime('2016-06-24'), pd.to_datetime('2016-07-22')]\n",
        "\n",
        "# Lookback window used for forecasting (usually 1-5x the forecast horizon)\n",
        "CONTEXT_WINDOW = 1 * POST_PERIOD_LENGTH"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "15c5KU1nv-Pd"
      },
      "source": [
        "# Other constants (no need to change these)\n",
        "\n",
        "SCENARIO = 'exchange-rates'\n",
        "TRAIN_CSV = 'train.csv'\n",
        "FORECAST_INPUT_CSV = 'forecast_input.csv'"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0z9SesuwAsqb"
      },
      "source": [
        "## Setup Google Cloud project"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tLpl54g9BSjk"
      },
      "source": [
        "# Authenticate to Google Cloud\n",
        "\n",
        "if 'google.colab' in sys.modules: \n",
        "    from google.colab import auth\n",
        "    auth.authenticate_user()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6cqbdePVAre9"
      },
      "source": [
        "# Initialize Vertex AI SDK\n",
        "\n",
        "from google.cloud import aiplatform\n",
        "\n",
        "aiplatform.init(project=PROJECT_ID, location=REGION)\n",
        "\n",
        "!gcloud config set project $PROJECT_ID"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4_zLQ5fs8JWj"
      },
      "source": [
        "## Read data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uGgRfGztZ0T3"
      },
      "source": [
        "We will import two related time-series from [FRED](https://fred.stlouisfed.org/) (Federal Reserve Economic Data). One is the history of the US Dollar : British Pound exchange rate, and the other is of the US Dollar : Euro.  "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nPI5TpG6N6V-"
      },
      "source": [
        "def read_fred_data(dataset, start_date, end_date):\n",
        "  url = f'https://fred.stlouisfed.org/graph/fredgraph.csv?id={dataset}&cosd={str(start_date)[:10]}&coed={str(end_date)[:10]}'\n",
        "  return pd.read_csv(url, index_col='DATE', parse_dates=True, na_values='.')\n",
        "\n",
        "df_DEXUSUK=read_fred_data('DEXUSUK', PRE_PERIOD[0], POST_PERIOD[1])\n",
        "df_DEXUSEU=read_fred_data('DEXUSEU', PRE_PERIOD[0], POST_PERIOD[1])\n",
        "\n",
        "# Merge each series into one dataframe\n",
        "df = pd.merge(left=df_DEXUSEU, left_on=df_DEXUSEU.index, right=df_DEXUSUK, right_on=df_DEXUSUK.index).rename(columns={'key_0':'DATE'}).set_index('DATE').dropna()\n",
        "\n",
        "df.head()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "L6sunUOE00gy"
      },
      "source": [
        "## Visualize data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VT5JaC6A2rO4"
      },
      "source": [
        "Let's now plot both time series: British Pound (`DEXUSUK`) and Euro (`DEXUSEU`).\n",
        "\n",
        "The shaded area is the four-week \"post period,\" which indicates a drop after the Brexit vote on June 23, 2016.\n",
        "\n",
        "We will later generate \"what would have been\" (counterfactual) by learning from the proxy Euro time series.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ycH9JaOFN5-q"
      },
      "source": [
        "fig, ax = plt.subplots()\n",
        "fig.set_figwidth(20)\n",
        "ax.plot(df[POST_PERIOD[1] - timedelta(days=730):POST_PERIOD[1]])\n",
        "ax.axvspan(POST_PERIOD[0], POST_PERIOD[1], alpha=0.5, color='lightblue')\n",
        "\n",
        "plt.title('GBP/US and Euro/US Exchange Rates')\n",
        "plt.legend(df.columns)\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uGicQZMU3nGn"
      },
      "source": [
        "Calculate the correlation between time series, to validate that an appropriate proxy is being used."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nkhn6f_OsCmJ"
      },
      "source": [
        "corr, _ = stats.pearsonr(x=df.DEXUSEU[df.index <= PRE_PERIOD[1]], y=df.DEXUSUK[df.index <= PRE_PERIOD[1]])\n",
        "\n",
        "print(f'Pearson correlation between time series in the pre-period: {round(corr, 2)}')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4_ikcYyyA4yk"
      },
      "source": [
        "## Create and upload data files"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nUbymt0n3013"
      },
      "source": [
        "We will now create and upload CSV files in the format that Vertex AI expects for training and forecasting."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EZ7UvtaeA6HY"
      },
      "source": [
        "# Training data\n",
        "\n",
        "df_train = df[df.index <= PRE_PERIOD[1]].copy()\n",
        "\n",
        "# There is only one series, so we will set the series ID to 0 for all rows.\n",
        "df_train['ID'] = 0\n",
        "\n",
        "df_train.to_csv(TRAIN_CSV)\n",
        "\n",
        "!gsutil cp {TRAIN_CSV} {STORAGE_PATH}/{TRAIN_CSV}"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XRI8fpwIEdpU"
      },
      "source": [
        "# Forecast input: create dataframe\n",
        "\n",
        "# Start with treatment dataframe\n",
        "df_forecast_input = df.copy()\n",
        "\n",
        "df_forecast_input['ID'] = 0\n",
        "\n",
        "# Set sales to empty during treatment period (to be predicted)\n",
        "df_forecast_input.loc[df_forecast_input.index >= POST_PERIOD[0], 'DEXUSUK'] = np.NaN\n",
        "\n",
        "df_forecast_input"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5C5eAWAMGg4W"
      },
      "source": [
        "# Forecast input: upload\n",
        "\n",
        "df_forecast_input.to_csv(FORECAST_INPUT_CSV)\n",
        "\n",
        "!gsutil cp {FORECAST_INPUT_CSV} {STORAGE_PATH}/{FORECAST_INPUT_CSV}"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "z6baKrrkBZwU"
      },
      "source": [
        "# Create dataset\n",
        "\n",
        "\n",
        "ds = dataset = aiplatform.TimeSeriesDataset.create(\n",
        "    display_name=f'dataset-{SCENARIO}',\n",
        "    gcs_source=f'{STORAGE_PATH}/{TRAIN_CSV}',\n",
        ")\n",
        "\n",
        "ds.resource_name"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kgLjGsjoEvub"
      },
      "source": [
        "# Train forecasting job\n",
        "\n",
        "job = aiplatform.AutoMLForecastingTrainingJob(\n",
        "    display_name=f'model-{SCENARIO}',\n",
        "    column_transformations=[\n",
        "        {\"numeric\": {\"column_name\": \"DEXUSUK\"}},\n",
        "        {\"timestamp\": {\"column_name\": \"DATE\"}},\n",
        "    ],\n",
        "    optimization_objective='minimize-rmse'\n",
        ")\n",
        "\n",
        "# This may take several hours to run\n",
        "model = job.run(\n",
        "    dataset=ds,\n",
        "    target_column='DEXUSUK',\n",
        "    time_column='DATE',\n",
        "    time_series_identifier_column='ID',\n",
        "    unavailable_at_forecast_columns=['DEXUSUK'],\n",
        "    available_at_forecast_columns=['DEXUSEU', 'DATE'],\n",
        "    data_granularity_unit='day',\n",
        "    forecast_horizon=POST_PERIOD_LENGTH,\n",
        "    data_granularity_count=1,\n",
        "    context_window=CONTEXT_WINDOW,\n",
        "    budget_milli_node_hours=8000\n",
        ")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "T3iyx-z5iXYL"
      },
      "source": [
        "## Forecast with trained model"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NaN8DlFQicez"
      },
      "source": [
        "# Try to access trained model from previous step.\n",
        "# If one is not available (e.g. because the runtime has timed out), you can update YOUR-MODEL-ID below to use an already trained model.\n",
        "\n",
        "try:\n",
        "  model\n",
        "except:\n",
        "  model = aiplatform.Model('YOUR-MODEL-ID') # Example: 7354146250200646848"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sRznGhsKjqC_"
      },
      "source": [
        "# Perform a batch prediction\n",
        "\n",
        "batch_prediction_job = model.batch_predict(\n",
        "    job_display_name=f'prediction-{SCENARIO}',\n",
        "    gcs_source=f'{STORAGE_PATH}/{FORECAST_INPUT_CSV}',\n",
        "    gcs_destination_prefix=f'{STORAGE_PATH}',\n",
        "    instances_format='csv',\n",
        "    predictions_format='csv'\n",
        ")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "aRe4yjoJVs6M"
      },
      "source": [
        "# The job output provides the full path with the timestamp where the output CSV was generated.\n",
        "# You can also look up this by looking at the batch prediction job details in the UI.\n",
        "\n",
        "try:\n",
        "  batch_prediction_job\n",
        "except:\n",
        "  batch_prediction_job = aiplatform.BatchPredictionJob('YOUR-BATCH-PREDICTION-JOB-ID')  # Example: 7296593272910163968\n",
        "\n",
        "FORECAST_OUTPUT_FILE = f'{batch_prediction_job.output_info.gcs_output_directory}/predictions_1.csv'\n",
        "\n",
        "df_forecast_output = pd.read_csv(FORECAST_OUTPUT_FILE, index_col='DATE', parse_dates=True).sort_values('DATE')\n",
        "\n",
        "df_forecast_output.head()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "brAusDx1do4M"
      },
      "source": [
        "# We will now create a dataframe that contains the actual and predicted values for comparison.\n",
        "\n",
        "df_estimate = df[['DEXUSUK']]\n",
        "\n",
        "df_counterfactual = df_estimate.copy().rename(columns={'DEXUSUK':'predicted_DEXUSUK'})\n",
        "df_counterfactual.loc[df_counterfactual.index.isin(df_estimate.index), 'predicted_DEXUSUK'] = df_forecast_output['predicted_DEXUSUK']\n",
        "\n",
        "df_estimate = pd.merge(df_estimate, df_counterfactual, on='DATE')\n",
        "df_estimate"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wc1i-dMO4vKS"
      },
      "source": [
        "## Plot the forecast"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "to_48Hbgh_Zy"
      },
      "source": [
        "fig, ax = plt.subplots()\n",
        "fig.set_figwidth(20)\n",
        "ax.plot(df_estimate[df_estimate.index >= POST_PERIOD[0] - timedelta(days=3 * POST_PERIOD_LENGTH)])\n",
        "ax.fill_between(df_estimate.index, df_estimate.DEXUSUK, df_estimate.predicted_DEXUSUK, color='lightblue')\n",
        "\n",
        "plt.title('Comparison of actual and predicted US/UK exchange rate')\n",
        "plt.legend(df_estimate.columns)\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "F7J6fA0t49Fd"
      },
      "source": [
        "## Calculate summary statistics"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "i-Q6fQNbjzhB"
      },
      "source": [
        "df_post_period = df_estimate[df_estimate.index >= POST_PERIOD[0]]\n",
        "\n",
        "mean_actual = df_post_period.DEXUSUK.mean()\n",
        "mean_predicted = df_post_period.predicted_DEXUSUK.mean()\n",
        "\n",
        "print(f'Actual US/UK exchange rate:    {round(mean_actual, 3)}') \n",
        "print(f'Predicted US/UK exchange rate: {round(mean_predicted, 3)}')\n",
        "\n",
        "auc = (df_post_period.DEXUSUK - df_post_period.predicted_DEXUSUK).sum()\n",
        "ate = auc / len(df_post_period.index)\n",
        "\n",
        "print(f'Average treatment effect:     {round(ate, 3)}')\n",
        "print(f'Relative treatment effect:    {round(ate/mean_predicted * 100, 3)}%')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mouId1dZVS3d"
      },
      "source": [
        "### Estimate effect with tfcausalimpact"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DUkmrjxzgDTy"
      },
      "source": [
        "We've already gone through data preprocessing steps, so it should be straightforward to pass in inputs into the [tfcausalimpact library](https://github.com/WillianFuks/tfcausalimpact), which uses [TensorFlow Probability](https://www.tensorflow.org/probability). It provides prediction intervals, to enable hypothesis testing on the significance of the treatment effect. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "N3VPJxKzVVhd"
      },
      "source": [
        "ci = CausalImpact(df[['DEXUSUK', 'DEXUSEU']], PRE_PERIOD, POST_PERIOD)\n",
        "print(ci.summary())\n",
        "print(ci.summary(output='report'))\n",
        "ci.plot()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GBlLAqPDVaTe"
      },
      "source": [
        "## Conclusion"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "J5xmGZytpG67"
      },
      "source": [
        "In this notebook, we walked through how to estimate the effect of an intervention. We used 2 methods to generate the counterfactual, or synthetic control, that allows us to compare the outcomes with and without the treatment."
      ]
    }
  ]
}