!pip install google-cloud-bigquery google-cloud-bigquery-storage --upgrade

# Automatically restart kernel after installs
import IPython

app = IPython.Application.instance()
app.kernel.do_shutdown(True)

"""
Explanation: <table align="left">
    <td>
      <a href="https://console.cloud.google.com/ai-platform/notebooks/deploy-notebook?name=BigQuery%20ML%20-%20Retail%20Demand%20Forecasting&download_url=https%3A%2F%2Fgithub.com%2FGoogleCloudPlatform%2Fanalytics-componentized-patterns%2Fraw%2Fmaster%2Fretail%2Ftime-series%2Fbqml-demand-forecasting%2Fbqml_retail_demand_forecasting.ipynb&url=https%3A%2F%2Fgithub.com%2FGoogleCloudPlatform%2Fanalytics-componentized-patterns%2Ftree%2Fmaster%2Fretail%2Ftime-series%2Fbqml-demand-forecasting%2F">
        <img src="https://cloud.google.com/images/products/ai/ai-solutions-icon.svg" alt="AI Platform Notebooks">Run on AI Platform Notebooks</a>
    </td>
  <td>
    <a href="https://github.com/GoogleCloudPlatform/analytics-componentized-patterns/tree/master/retail/time-series/bqml-demand-forecasting/">
      <img src="https://cloud.google.com/ml-engine/images/github-logo-32px.png" alt="GitHub logo">
      View on GitHub
    </a>
  </td>
</table>

Building a demand forecasting model by using BigQuery ML
This notebook shows you how to train, deploy, and evaluate a time series model by using BigQuery ML. It provides an end-to-end solution for forecasting multiple products. Using the public Iowa Liquor Sales data dataset, you use a single SQL query to create five time series models, where each model forecasts the retail sales of a single liquor product. 
By the end of this notebook, you will know how to:
* Pre-process time series data into the correct format needed to create the model.
* Train the time series model in BigQuery ML.
* Evaluate the model.
* Make predictions about future demand using the model.
* Create a dashboard to visualize the forecasted demand using Data Studio.
Setup
Install the required packages, configure the environment variables, and create a BigQuery dataset.
PIP Install Packages and dependencies
End of explanation
"""


PROJECT_ID = "yourProject"  # Change to your project.
REGION = "US"

"""
Explanation: Configure GCP environment settings
Update the following variables to reflect the values for your GCP environment:

PROJECT_ID: The ID of the Google Cloud project you are using to implement this solution.
REGION: The region to use for the BigQuery dataset.
End of explanation
"""


from datetime import datetime, timedelta

import pandas as pd
from google.cloud import bigquery

pd.set_option("display.float_format", lambda x: "%.3f" % x)

"""
Explanation: Import libraries and define constants
End of explanation
"""


import matplotlib.pyplot as plt
import pandas as pd


def plot_historical_and_forecast(
    input_timeseries,
    timestamp_col_name,
    data_col_name,
    forecast_output=None,
    actual=None,
    title=None,
    plotstartdate=None,
):

    if plotstartdate:
        input_timeseries = input_timeseries[
            input_timeseries[timestamp_col_name] >= pd.to_datetime(plotstartdate)
        ]
    input_timeseries = input_timeseries.sort_values(timestamp_col_name)

    # Plot the input historical data
    plt.figure(figsize=(20, 6))
    plt.plot(
        input_timeseries[timestamp_col_name],
        input_timeseries[data_col_name],
        label="Historical",
    )
    plt.xlabel(timestamp_col_name)
    plt.ylabel(data_col_name)

    if forecast_output is not None:
        forecast_output = forecast_output.sort_values("forecast_timestamp")
        forecast_output["forecast_timestamp"] = pd.to_datetime(
            forecast_output["forecast_timestamp"]
        )
        x_data = forecast_output["forecast_timestamp"]
        y_data = forecast_output["forecast_value"]
        confidence_level = forecast_output["confidence_level"].iloc[0] * 100
        low_CI = forecast_output["confidence_interval_lower_bound"]
        upper_CI = forecast_output["confidence_interval_upper_bound"]
        # Plot the forecast data
        plt.plot(x_data, y_data, alpha=1, label="Forecast", linestyle="--")
        # Shade the confidence interval
        plt.fill_between(
            x_data,
            low_CI,
            upper_CI,
            color="#539caf",
            alpha=0.4,
            label=f"{confidence_level} confidence interval",
        )

    # Plot actual data
    if actual is not None:
        actual = actual.sort_values(timestamp_col_name)
        plt.plot(
            actual[timestamp_col_name],
            actual[data_col_name],
            label="Actual",
            linestyle="--",
        )

    # Display title, legend
    plt.title(f"{title}", fontsize=16)
    plt.legend(loc="upper center", prop={"size": 16})

"""
Explanation: Import libraries for data manipulation and plotting
End of explanation
"""


!bq mk --location=$REGION --dataset $PROJECT_ID:bqmlforecast

"""
Explanation: Create a BigQuery dataset
Create a dataset in your project called bqmlforecast:
End of explanation
"""


%%bigquery --project $PROJECT_ID

SELECT 
    invoice_and_item_number,
    date,
    store_number,
    item_description,
    bottles_sold,
    sale_dollars
FROM
  `bigquery-public-data.iowa_liquor_sales.sales` 
LIMIT 
  5

"""
Explanation: Prepare the training data
You train the time series models on a dataset containing transactional data. Each row represents a transaction on a single product, as identified by the item_description value, and contains details such as the number of bottles sold and the sales amount in dollars. In following steps, you use the number of bottles sold value to forecast product demand.
Note: Jupyter runs cells starting with %%bigquery as SQL queries
End of explanation
"""


# %%bigquery --project $PROJECT_ID

# CREATE OR REPLACE VIEW bqmlforecast.training_data AS (
# SELECT
#   date,
#   item_name,
#   total_amount_sold
# FROM
#   `[YOUR_PROJECT].[YOUR_DATASET].[YOUR_SOURCE_TABLE]`
# );

"""
Explanation: [Optional] Transform your data to the expected schema
If you use the sample dataset, you can skip this step.
If you want to use your own data, import your sales data to BigQuery, and then adapt the following SQL query to create a view with the appropriate schema. Update [YOUR_PROJECT].[YOUR_DATASET].[YOUR_SOURCE_TABLE] to reflect the project, dataset, and table containing your sales data.
End of explanation
"""


ARIMA_PARAMS = {
    "TRAININGDATA_STARTDATE": "2016-01-01",
    "TRAININGDATA_ENDDATE": "2017-06-01",
}
ARIMA_PARAMS

"""
Explanation: Set the start and end dates for the training data
You can adjust the TRAININGDATA_STARTDATE and TRAININGDATA_ENDDATE parameters to specify the start/end dates of your training data:
End of explanation
"""


%%bigquery --params $ARIMA_PARAMS  --project $PROJECT_ID 

CREATE OR REPLACE TABLE bqmlforecast.training_data AS (
    WITH topsellingitems AS(
         SELECT 
            item_description,
            count(item_description) cnt_transactions
        FROM
            `bigquery-public-data.iowa_liquor_sales.sales` 
        GROUP BY 
            item_description
        ORDER BY cnt_transactions DESC
        LIMIT 5 #Top N
    )
    SELECT 
        date,
        item_description AS item_name,
        SUM(bottles_sold) AS total_amount_sold
    FROM
        `bigquery-public-data.iowa_liquor_sales.sales` 
    GROUP BY
        date, item_name
    HAVING 
        date BETWEEN @TRAININGDATA_STARTDATE AND @TRAININGDATA_ENDDATE
        AND item_description IN (SELECT item_description FROM topsellingitems)
    );

SELECT 
    date,
    item_name,
    total_amount_sold
FROM 
    bqmlforecast.training_data 
ORDER BY item_name, date
LIMIT 10

"""
Explanation: Write the training data to a table
If you review the training data, you can observe that some days have no transactions for a given product. To keep you from having to do additional pre-processing, BigQuery ML automatically handles the following situations:

Missing values: Values are imputed using local linear interpolation.
Duplicated timestamps: Values are averaged across duplicated timestamps.
Spike and dip anomalies: Values are standardized using local z-scores.
End of explanation
"""


%%bigquery dfhistorical --project $PROJECT_ID 

SELECT 
    * 
FROM 
    bqmlforecast.training_data

"""
Explanation: Plot the sales histories of the target liquor products
Save the training data to the dfhistorical Pandas dataframe:
End of explanation
"""


itemslist = list(dfhistorical.item_name.unique())

for item in itemslist:

    datah = dfhistorical[dfhistorical.item_name == item]
    plot_historical_and_forecast(
        input_timeseries=datah,
        timestamp_col_name="date",
        data_col_name="total_amount_sold",
        forecast_output=None,
        actual=None,
        title=item,
    )

"""
Explanation: Using the training data, plot the sales histories of the target liquor products:
End of explanation
"""


%%bigquery --project $PROJECT_ID 

CREATE OR REPLACE MODEL bqmlforecast.arima_model

OPTIONS(
  MODEL_TYPE='ARIMA',
  TIME_SERIES_TIMESTAMP_COL='date', 
  TIME_SERIES_DATA_COL='total_amount_sold',
  TIME_SERIES_ID_COL='item_name',
  HOLIDAY_REGION='US'
) AS

SELECT 
    date,
    item_name,
    total_amount_sold
FROM
  bqmlforecast.training_data

"""
Explanation: Train the model
Because you are training the model on multiple products in a single model creation statement, you must specify the item_name column for the TIME_SERIES_ID_COL parameter. If you were only forecasting a single item, then you would not need to specify TIME_SERIES_ID_COL. For more information, see the BigQuery ML time series model creation documentation.
You can also account for holiday effects when doing time series modeling in BigQuery ML. By default, holiday effects modeling is disabled. But since this data is from the United States, and the data includes a minimum one year of daily data, you can also specify an optional HOLIDAY_REGION. With holiday effects enabled, spikes and dips that appear during holidays will no longer be treated as anomalies. A full list of the holiday regions can be found in the HOLIDAY_REGION documentation.
End of explanation
"""


%%bigquery --project $PROJECT_ID 

SELECT
  *
FROM
  ML.EVALUATE(MODEL bqmlforecast.arima_model)

"""
Explanation: Evaluate the model
Use the ML.EVALUATE function to see the evaluation metrics of all the created models. 
The non_seasonal_{p,d,q} and has_drift columns define the time series model. The log_likelihood, AIC, and variance columns are relevant to the model fitting process. The fitting process determines the best model by using the auto.ARIMA algorithm, one for each time series.
End of explanation
"""


%%bigquery dfforecast --project $PROJECT_ID 

DECLARE HORIZON STRING DEFAULT "30"; #number of values to forecast
DECLARE CONFIDENCE_LEVEL STRING DEFAULT "0.90";

EXECUTE IMMEDIATE format("""
    SELECT
      *
    FROM 
      ML.FORECAST(MODEL bqmlforecast.arima_model, 
                  STRUCT(%s AS horizon, 
                         %s AS confidence_level)
                 )
    """,HORIZON,CONFIDENCE_LEVEL)

dfforecast.head()

"""
Explanation: You can see that five models were trained, one for each of the products in the item_name column. Each model has its own p,d,q hyperparameters, and the detected seasonality for these five models is WEEKLY.
Make predictions using the model
Make predictions using ML.FORECAST, which forecasts the next n values, as specified in the horizon parameter. You can also optionally change the confidence_level parameter to change the percentage of future values that fall in the prediction interval. You save the prediction data to the dfforecast dataframe so that you can plot it in a following step.
End of explanation
"""


print(f"Number of rows: {dfforecast.shape[0]}")

"""
Explanation: Since horizon is set to 30, the result is 30 x (number of items), with one row per forecasted value:
End of explanation
"""


%%bigquery --project $PROJECT_ID 

SELECT
  *
FROM 
  ML.ARIMA_COEFFICIENTS(MODEL bqmlforecast.arima_model)

"""
Explanation: Inspect the model coefficients
You can view the coefficients of each of the time series models using ML.ARIMA_COEFFICIENTS.
For each of the models, the ar_coefficients value shows the model coefficients of the autoregressive (AR) part of the time series model. Similarly, the ma_coefficients value shows the model coefficients of the moving-average (MA) part of the model. Both of these values are arrays, and their lengths are equal to non_seasonal_p and non_seasonal_q, respectively. The intercept_or_drift value is the constant term in the time series model.
End of explanation
"""


itemslist = list(dfhistorical.item_name.unique())

for item in itemslist:
    datah = dfhistorical[dfhistorical.item_name == item]
    dataf = dfforecast[dfforecast.item_name == item]

    plot_historical_and_forecast(
        input_timeseries=datah,
        timestamp_col_name="date",
        data_col_name="total_amount_sold",
        forecast_output=dataf,
        actual=None,
        title=item,
        plotstartdate="2017-01-01",
    )

"""
Explanation: Plot the forecasted predictions against the historical data
Plot the forecasted predictions, using the dfhistorical dataframe that contains the historical data that you used for training and the dfforecast dataframe that contains the prediction data:
End of explanation
"""


%%bigquery dfactual --params $ARIMA_PARAMS --project $PROJECT_ID 

DECLARE HORIZON STRING DEFAULT "30"; #number of values to forecast

SELECT 
    date,
    item_description AS item_name,
    SUM(bottles_sold) AS total_amount_sold
FROM
    `bigquery-public-data.iowa_liquor_sales.sales` 
GROUP BY
    date, item_name
HAVING 
    date BETWEEN DATE_ADD(@TRAININGDATA_ENDDATE, 
                              INTERVAL 1 DAY) 
            AND DATE_ADD(@TRAININGDATA_ENDDATE, 
                             INTERVAL 1+CAST(HORIZON AS INT64) DAY) 
ORDER BY
    date;

dfactual.head()

"""
Explanation: Plot the forecasted predictions against the actual data
Save the actual sales data to the dfactual Pandas dataframe:
End of explanation
"""


itemslist = list(dfhistorical.item_name.unique())

for item in itemslist:
    datah = dfhistorical[dfhistorical.item_name == item].sort_values("date")
    dataf = dfforecast[dfforecast.item_name == item].sort_values(["forecast_timestamp"])
    dataa = dfactual[dfactual.item_name == item].sort_values("date")
    plot_historical_and_forecast(
        input_timeseries=datah,
        timestamp_col_name="date",
        data_col_name="total_amount_sold",
        forecast_output=dataf,
        actual=dataa,
        title=item,
        plotstartdate="2017-01-01",
    )

"""
Explanation: Plot the forecasted predictions against the training data and actual values, using the dfhistorical dataframe that contains the historical data, the dfforecast dataframe that contains the prediction data, and the dfactual dataframe that contains the actual sales data:
End of explanation
"""


%%bigquery --params $ARIMA_PARAMS  --project $PROJECT_ID 

CREATE OR REPLACE VIEW bqmlforecast.outputdata_datastudio AS (
  SELECT
    date AS timestamp,
    item_name,
    total_amount_sold AS history_value,
    NULL AS forecast_value,
    NULL AS prediction_interval_lower_bound,
    NULL AS prediction_interval_upper_bound
  FROM
    bqmlforecast.training_data
  UNION ALL
  SELECT
    EXTRACT(DATE
    FROM
      forecast_timestamp) AS timestamp,
    item_name,
    NULL AS history_value,
    forecast_value,
    prediction_interval_lower_bound,
    prediction_interval_upper_bound
  FROM
    ML.FORECAST(MODEL bqmlforecast.arima_model,
      STRUCT(30 AS horizon, 0.9 AS confidence_level)) 
  ORDER BY timestamp
  )

"""
Explanation: Create a dashboard with Data Studio
Follow the steps in this section to create an interactive, shareable dashboard of the forecasted data by using Data Studio. 
Create a view containing the data for the dashboard
Create a view that concatenates the historical and forecasted data. The SQL before the UNION ALL clause selects the historical data. The SQL after the UNION ALL clause uses ML.FORECAST to generate the forecasted value and the prediction interval. The query uses different fields for history_value and forecasted_value so that you can plot them in different colors in a following step.
End of explanation
"""


