{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# plotting libraries\nimport matplotlib\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom statsmodels.tsa.api import ExponentialSmoothing,SimpleExpSmoothing, Holt\nfrom matplotlib.dates import (\n        MonthLocator,\n        num2date,\n        AutoDateLocator,\n        AutoDateFormatter,\n)\nimport gc # garbage collector\n\n# stats models\nimport statsmodels.api as sm\nfrom fbprophet import Prophet\n\n# time libraries\nimport datetime\n\n# warning libraries for debugging\nimport warnings\n\n# deal with date in x-axis of plots\nfrom pandas.plotting import register_matplotlib_converters\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('../input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"code to create time bar to run functions"},{"metadata":{"trusted":true},"cell_type":"code","source":"import time, sys\nfrom IPython.display import clear_output\n\ndef update_progress(progress):\n    bar_length = 20\n    if isinstance(progress, int):\n        progress = float(progress)\n    if not isinstance(progress, float):\n        progress = 0\n    if progress < 0:\n        progress = 0\n    if progress >= 1:\n        progress = 1\n\n    block = int(round(bar_length * progress))\n\n    clear_output(wait = True)\n    text = \"Progress: [{0}] {1:.1f}%\".format( \"#\" * block + \"-\" * (bar_length - block), progress * 100)\n    print(text)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","collapsed":true,"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":false},"cell_type":"markdown","source":"# Introduction\nI am still a novice when it comes to time series so I want to start simple. In this kernel we are interested in predicting meter values based on our time series data. We are going to ignore all the other features. Let's jump in..."},{"metadata":{},"cell_type":"markdown","source":"# Load Data\nimport the training dataframe"},{"metadata":{"trusted":true},"cell_type":"code","source":"%time train_df = pd.read_csv('../input/ashrae-energy-prediction/train.csv')\n%time test_df = pd.read_csv('../input/ashrae-energy-prediction/test.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Set Column Datatypes"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Saving some memory\nd_types = {'building_id': np.int16,\n          'meter': np.int8}\n\nfor feature in d_types:\n    train_df[feature] = train_df[feature].astype(d_types[feature])\n    \n    \ntrain_df[\"timestamp\"] = pd.to_datetime(train_df[\"timestamp\"], infer_datetime_format=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Add log of meter values"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df[\"log_meter_reading\"]=np.log(train_df[\"meter_reading\"]+.00001)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Evaluation Metric\nThe evaluation metric for this competition is the root mean squared logarithmic error. Below I created a method that can calculate this value."},{"metadata":{"trusted":true},"cell_type":"code","source":"def rmsle(pred_series,true_series):\n    sum_series = (np.log(pred_series+1) - \\\n        np.log(true_series+1))**2\n    return np.sqrt(np.sum(sum_series))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Splitting the Training Dataset\nWe need to split the original training data into a training and validation set. I decided to splits the training set into the first 9 months and the validation set into the last 3 months."},{"metadata":{"trusted":true},"cell_type":"code","source":"start_validation='2016-12-15'\ntrain = train_df.loc[train_df[\"timestamp\"]<start_validation,:]\nvalid = train_df.loc[train_df[\"timestamp\"]>=start_validation,:]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Some algorithms need the training data to be indexed by the time stamps."},{"metadata":{"trusted":true},"cell_type":"code","source":"# for the training data I want to reformat\n# the dataframe so that the timestamp is the \n# index\nprint(\"reformat training data frame...\")\ndef trainDF2timeDF(training_df):\n    timeValue_df =  training_df.copy()\n    timeValue_df = timeValue_df.set_index(\"timestamp\")\n    warnings.simplefilter(\"ignore\")\n    timeValue_df.index = pd.to_datetime(timeValue_df.index.values)\n    return(timeValue_df)\n\ntimeIndexed_train = trainDF2timeDF(train)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Naive Approach: Just take the average values from the training data\nFor this approach I use the average of the training data to predict the validation data. I will use this method as my baseline for more complicated modeling."},{"metadata":{"trusted":true},"cell_type":"code","source":"# create new data frame for this model\nvalid_avgVal_df = valid.copy()\n# rename timestamp to signify the current meter reading time\nvalid_avgVal_df = valid_avgVal_df.rename(\n    columns={\"timestamp\": \"now\", \n             \"meter_reading\": \"cur_meter_reading\",\n            \"log_meter_reading\":\"cur_log_meter_reading\"})\n\n# This model splits the data based on \n# building ID and model type\nnbuildings=len(valid[\"building_id\"].unique())\nprint(\"number of buildings: \"+ str(nbuildings))\nx=0\nfor b_id in list(valid[\"building_id\"].unique()):\n    update_progress(x / nbuildings)\n    x+=1\n    for meter_t in list(\n        valid_avgVal_df.loc[valid_avgVal_df[\"building_id\"]==b_id,\"meter\"].unique()):\n        if(not ((b_id in train[\"building_id\"]) and\n           (meter_t in train.loc[train[\"building_id\"]==b_id,\"meter\"].values))):\n            print(\"missing!\")\n            print(b_id)\n            print(meter_t)\n            # if there is no meter reading for a specific\n            # building ID then I'll just set the reading\n            # to the average value of that meter given\n            # all of the building IDs.\n            valid.loc[((valid_avgVal_df[\"building_id\"]==b_id) &\n                valid_avgVal_df[\"meter\"]==meter_t),\"pred_meter_reading\"] = \\\n                train.loc[train[\"meter\"]==meter_t,\"meter_reading\"].mean()\n        else:\n            # calculate the average meter_reading values\n            # for each meter given the building id\n            valid_avgVal_df.loc[((valid_avgVal_df[\"building_id\"]==b_id) &\n                valid_avgVal_df[\"meter\"]==meter_t),\"pred_meter_reading\"] = \\\n                train.loc[(\n                (train[\"building_id\"]==b_id) &\n                (train[\"meter\"]==meter_t)),\"meter_reading\"].mean()\nupdate_progress(1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can visual our predictions for each building and meter as such"},{"metadata":{"trusted":true},"cell_type":"code","source":"b_i=1\nm_t=0\ntrain_bidX_meterY = train.loc[(\n    (train[\"building_id\"]==b_i) &\n    (train[\"meter\"]==m_t)),:].copy()\nvalid_bidX_meterY = valid.loc[(\n    (valid[\"building_id\"]==b_i) &\n    (valid[\"meter\"]==m_t)),:].copy()\npred_bidX_meterY = valid_avgVal_df.loc[(\n    (valid_avgVal_df[\"building_id\"]==b_i) &\n    (valid_avgVal_df[\"meter\"]==m_t)),:].copy()\n\nplt.figure(figsize =(15,8))\nplt.plot(train_bidX_meterY['meter_reading'], label = 'Train')\nplt.plot(valid_bidX_meterY['meter_reading'], label = 'Validation')\nplt.plot(pred_bidX_meterY['pred_meter_reading'], label = 'Simple Exponential Smoothing')\nplt.legend(loc = 'best')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now let's evaluate our model."},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Naive Approach - RMSLE value:\")\nprint(rmsle(valid_avgVal_df[\"pred_meter_reading\"],\n           valid_avgVal_df[\"cur_meter_reading\"]))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This is a really high value but it provides a good baseline moving forward."},{"metadata":{},"cell_type":"markdown","source":"Since each meter seems to have its own patter I am also interested in the RMSLE value for each meter."},{"metadata":{"trusted":true},"cell_type":"code","source":"avgVal_rmsle_list=[]\nfor meter_t in list(valid_avgVal_df[\"meter\"].unique()):\n        sub_valid_avgVal_df = valid_avgVal_df.loc[(\n            valid_avgVal_df[\"meter\"]==meter_t),:].copy()\n        sub_rmsle = rmsle(sub_valid_avgVal_df[\"pred_meter_reading\"],\n           sub_valid_avgVal_df[\"cur_meter_reading\"])\n        sub_rmsle_df = pd.DataFrame({\"meter\":[meter_t],\n                                   \"rmsle\":[sub_rmsle]})\n        avgVal_rmsle_list.append(sub_rmsle_df)\navgVal_rmsle_df = pd.concat(avgVal_rmsle_list)\navgVal_rmsle_df  ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Simple Exponential Smoothing\nIn our baseline model, we took the average of past meter values to predict the future meter values. However, instead of weighing each past meter value equally we can assume that the most recent reading should probably way higher than readings from the distant past. **Simple Exponential Smoothing** uses weighted averages to give the largest weights to the most recent observations and the smallest weights to the oldest observations. The weights or \"smoothing\" parameter have been labeled as $\\alpha$. So our forcast now looks like:\n$\\hat{y}_{T+1|T}=\\alpha y_{T} + \\alpha(1-\\alpha)y_{T-1} + \\alpha(1-\\alpha)^{2}y_{T-2}+...,$"},{"metadata":{},"cell_type":"markdown","source":"Let's try this method. We can use the auto optimization to automatically find an optimized $\\alpha$ value for us."},{"metadata":{"trusted":true},"cell_type":"code","source":"valid_expSmooth = valid.copy()\n# rename timestamp to signify the current meter reading time\nvalid_expSmooth = valid_expSmooth.rename(\n    columns={\"timestamp\": \"now\", \n             \"meter_reading\": \"cur_meter_reading\",\n            \"log_meter_reading\":\"cur_log_meter_reading\"})\n\n# for the training data I want to reformat\n# the dataframe so that the timestamp is the \n# index\nprint(\"reformat training data frame...\")\ndef trainDF2timeDF(training_df):\n    timeValue_df =  train.copy()\n    timeValue_df = timeValue_df.set_index(\"timestamp\")\n    warnings.simplefilter(\"ignore\")\n    timeValue_df.index = pd.to_datetime(timeValue_df.index.values)\n    return(timeValue_df)\n\ntimeIndexed_train = trainDF2timeDF(train)\n\n# This model splits the data based on \n# building ID and model type\nnbuildings=len(valid[\"building_id\"].unique())\nprint(\"number of buildings: \"+ str(nbuildings))\nx=0\nfor b_id in list(valid[\"building_id\"].unique()):\n    update_progress(x / nbuildings)\n    x+=1\n    for meter_t in list(\n        valid_expSmooth.loc[valid_expSmooth[\"building_id\"]==b_id,\"meter\"].unique()):\n        if(not ((b_id in train[\"building_id\"]) and\n           (meter_t in train.loc[train[\"building_id\"]==b_id,\"meter\"].values))):\n            print(\"missing!\")\n            print(b_id)\n            print(meter_t)\n            # if there is no meter reading for a specific\n            # building ID then I'll just train\n            # independent of the building ID\n            sub_timeTrain_df = timeIndexed_train.loc[(\n                timeIndexed_train[\"meter\"]==meter_t),\"meter_reading\"].copy()\n            numValid = len(valid_expSmooth.loc[(\n                (valid_expSmooth[\"building_id\"]==b_id) &\n                (valid_expSmooth[\"meter\"]==meter_t)),:])\n            fit_simExpSmooth = SimpleExpSmoothing(sub_timeTrain_df).fit()\n            # forecast the meter_readings\n            valid_expSmooth.loc[(\n                (valid_expSmooth[\"building_id\"]==b_id) &\n                (valid_expSmooth[\"meter\"]==meter_t)),\"pred_meter_reading\"] = \\\n                fit_simExpSmooth.forecast(numValid).values\n            # collect the alpha level used\n            valid_expSmooth.loc[(\n                (valid_expSmooth[\"building_id\"]==b_id) &\n                (valid_expSmooth[\"meter\"]==meter_t)),\"alpha\"] = \\\n                fit_simExpSmooth.model.params['smoothing_level']\n        else:\n            # fit the model to the meter values of\n            # this building type\n            sub_timeTrain_df = timeIndexed_train.loc[(\n                (timeIndexed_train[\"building_id\"]==b_id) &\n                (timeIndexed_train[\"meter\"]==meter_t)),\"meter_reading\"].copy()\n            numValid = len(valid_expSmooth.loc[(\n                (valid_expSmooth[\"building_id\"]==b_id) &\n                (valid_expSmooth[\"meter\"]==meter_t)),:])\n            fit_simExpSmooth = SimpleExpSmoothing(sub_timeTrain_df).fit()\n            # forecast the meter_readings\n            valid_expSmooth.loc[(\n                (valid_expSmooth[\"building_id\"]==b_id) &\n                (valid_expSmooth[\"meter\"]==meter_t)),\"pred_meter_reading\"] = \\\n                fit_simExpSmooth.forecast(numValid).values\n            # collect the alpha level used\n            valid_expSmooth.loc[(\n                (valid_expSmooth[\"building_id\"]==b_id) &\n                (valid_expSmooth[\"meter\"]==meter_t)),\"alpha\"] = \\\n                fit_simExpSmooth.model.params['smoothing_level']\nupdate_progress(1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's see how it looks compared to the true values in building 0 meter 0 "},{"metadata":{"trusted":true},"cell_type":"code","source":"b_i=0\nm_t=0\ntrain_bidX_meterY = train.loc[(\n    (train[\"building_id\"]==b_i) &\n    (train[\"meter\"]==m_t)),:].copy()\nvalid_bidX_meterY = valid.loc[(\n    (valid[\"building_id\"]==b_i) &\n    (valid[\"meter\"]==m_t)),:].copy()\npred_bidX_meterY = valid_expSmooth.loc[(\n    (valid_expSmooth[\"building_id\"]==b_i) &\n    (valid_expSmooth[\"meter\"]==m_t)),:].copy()\n\nplt.figure(figsize =(15,8))\nplt.plot(train_bidX_meterY['meter_reading'], label = 'Train')\nplt.plot(valid_bidX_meterY['meter_reading'], label = 'Validation')\nplt.plot(pred_bidX_meterY['pred_meter_reading'], label = 'Simple Exponential Smoothing')\nplt.legend(loc = 'best')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Cool, now let me see some numbers..."},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Simple Exponential Smoothing - RMSLE value:\")\nprint(rmsle(valid_expSmooth[\"pred_meter_reading\"],\n           valid_expSmooth[\"cur_meter_reading\"]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"expSmooth_rmsle_list=[]\nfor meter_t in list(valid_expSmooth[\"meter\"].unique()):\n        sub_valid_expSmooth_df = valid_expSmooth.loc[(\n            valid_expSmooth[\"meter\"]==meter_t),:].copy()\n        sub_rmsle = rmsle(sub_valid_expSmooth_df[\"pred_meter_reading\"],\n           sub_valid_expSmooth_df[\"cur_meter_reading\"])\n        sub_rmsle_df = pd.DataFrame({\"meter\":[meter_t],\n                                   \"rmsle\":[sub_rmsle]})\n        expSmooth_rmsle_list.append(sub_rmsle_df)\nexpSmooth_rmsle_df = pd.concat(expSmooth_rmsle_list)\nexpSmooth_rmsle_df  ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Holt Model\nSo the last model did a lot better since it weighed the more recent meter values higher than the past ones instead of weighing them all equally. The problem is that as we get further from our nearest point in the training set, our meter readings should not stay the same. **Holt's model** accounts for trend. For example, it accounts for the fact that the meter readings may be going up/down/same over time. Holt's model is useful if there is no seasonality in the time series."},{"metadata":{"trusted":true},"cell_type":"code","source":"# create new data frame for this model\nvalid_holt = valid.copy()\n# rename timestamp to signify the current meter reading time\nvalid_holt = valid_holt.rename(\n    columns={\"timestamp\": \"now\", \n             \"meter_reading\": \"cur_meter_reading\",\n            \"log_meter_reading\":\"cur_log_meter_reading\"})\n\n\n# This model splits the data based on \n# building ID and model type\nnbuildings=len(valid[\"building_id\"].unique())\nprint(\"number of buildings: \"+ str(nbuildings))\nx=0\nfor b_id in list(valid[\"building_id\"].unique()):\n    update_progress(x / nbuildings)\n    x+=1\n    for meter_t in list(\n        valid_holt.loc[valid_holt[\"building_id\"]==b_id,\"meter\"].unique()):\n        if(not ((b_id in train[\"building_id\"]) and\n           (meter_t in train.loc[train[\"building_id\"]==b_id,\"meter\"].values))):\n            print(\"missing!\")\n            print(b_id)\n            print(meter_t)\n            # if there is no meter reading for a specific\n            # building ID then I'll just train\n            # independent of the building ID\n            sub_timeTrain_df = timeIndexed_train.loc[(\n                timeIndexed_train[\"meter\"]==meter_t),\"meter_reading\"].copy()\n            numValid = len(valid_holt.loc[(\n                (valid_holt[\"building_id\"]==b_id) &\n                (valid_holt[\"meter\"]==meter_t)),:])\n            fit_holt = Holt(\n                sub_timeTrain_df).fit(optimized=True)\n            # forecast the meter_readings\n            valid_holt.loc[(\n                (valid_holt[\"building_id\"]==b_id) &\n                (valid_holt[\"meter\"]==meter_t)),\"pred_meter_reading\"] = \\\n                fit_holt.forecast(numValid).values\n            # collect the alpha level used\n            valid_holt.loc[(\n                (valid_holt[\"building_id\"]==b_id) &\n                (valid_holt[\"meter\"]==meter_t)),\"alpha\"] = \\\n                fit_holt.model.params['smoothing_level']\n        else:\n            # fit the model to the meter values of\n            # this building type\n            sub_timeTrain_df = timeIndexed_train.loc[(\n                (timeIndexed_train[\"building_id\"]==b_id) &\n                (timeIndexed_train[\"meter\"]==meter_t)),\"meter_reading\"].copy()\n            numValid = len(valid_holt.loc[(\n                (valid_holt[\"building_id\"]==b_id) &\n                (valid_holt[\"meter\"]==meter_t)),:])\n            fit_holt = Holt(\n                sub_timeTrain_df).fit(optimized=True)\n            # forecast the meter_readings\n            valid_holt.loc[(\n                (valid_holt[\"building_id\"]==b_id) &\n                (valid_holt[\"meter\"]==meter_t)),\"pred_meter_reading\"] = \\\n                fit_holt.forecast(numValid).values\n            # collect the alpha level used\n            valid_holt.loc[(\n                (valid_holt[\"building_id\"]==b_id) &\n                (valid_holt[\"meter\"]==meter_t)),\"alpha\"] = \\\n                fit_holt.model.params['smoothing_level']\nupdate_progress(1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's see how it looks in each type of meter. I randomly chose building ID's with each type of meter"},{"metadata":{"trusted":true},"cell_type":"code","source":"b_i=0\nm_t=0\ntrain_bidX_meterY = train.loc[(\n    (train[\"building_id\"]==b_i) &\n    (train[\"meter\"]==m_t)),:].copy()\nvalid_bidX_meterY = valid.loc[(\n    (valid[\"building_id\"]==b_i) &\n    (valid[\"meter\"]==m_t)),:].copy()\npred_bidX_meterY = valid_holt.loc[(\n    (valid_holt[\"building_id\"]==b_i) &\n    (valid_holt[\"meter\"]==m_t)),:].copy()\n\nplt.figure(figsize =(15,8))\nplt.plot(train_bidX_meterY['meter_reading'], label = 'Train')\nplt.plot(valid_bidX_meterY['meter_reading'], label = 'Validation')\nplt.plot(pred_bidX_meterY['pred_meter_reading'], label = 'Holt Model')\nplt.legend(loc = 'best')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"b_i=161\nm_t=1\ntrain_bidX_meterY = train.loc[(\n    (train[\"building_id\"]==b_i) &\n    (train[\"meter\"]==m_t)),:].copy()\nvalid_bidX_meterY = valid.loc[(\n    (valid[\"building_id\"]==b_i) &\n    (valid[\"meter\"]==m_t)),:].copy()\npred_bidX_meterY = valid_holt.loc[(\n    (valid_holt[\"building_id\"]==b_i) &\n    (valid_holt[\"meter\"]==m_t)),:].copy()\n\nplt.figure(figsize =(15,8))\nplt.plot(train_bidX_meterY['meter_reading'], label = 'Train')\nplt.plot(valid_bidX_meterY['meter_reading'], label = 'Validation')\nplt.plot(pred_bidX_meterY['pred_meter_reading'], label = 'Holt Model')\nplt.legend(loc = 'best')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"b_i=745\nm_t=2\ntrain_bidX_meterY = train.loc[(\n    (train[\"building_id\"]==b_i) &\n    (train[\"meter\"]==m_t)),:].copy()\nvalid_bidX_meterY = valid.loc[(\n    (valid[\"building_id\"]==b_i) &\n    (valid[\"meter\"]==m_t)),:].copy()\npred_bidX_meterY = valid_holt.loc[(\n    (valid_holt[\"building_id\"]==b_i) &\n    (valid_holt[\"meter\"]==m_t)),:].copy()\n\nplt.figure(figsize =(15,8))\nplt.plot(train_bidX_meterY['meter_reading'], label = 'Train')\nplt.plot(valid_bidX_meterY['meter_reading'], label = 'Validation')\nplt.plot(pred_bidX_meterY['pred_meter_reading'], label = 'Holt Model')\nplt.legend(loc = 'best')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"b_i=106\nm_t=3\ntrain_bidX_meterY = train.loc[(\n    (train[\"building_id\"]==b_i) &\n    (train[\"meter\"]==m_t)),:].copy()\nvalid_bidX_meterY = valid.loc[(\n    (valid[\"building_id\"]==b_i) &\n    (valid[\"meter\"]==m_t)),:].copy()\npred_bidX_meterY = valid_holt.loc[(\n    (valid_holt[\"building_id\"]==b_i) &\n    (valid_holt[\"meter\"]==m_t)),:].copy()\n\nplt.figure(figsize =(15,8))\nplt.plot(train_bidX_meterY['meter_reading'], label = 'Train')\nplt.plot(valid_bidX_meterY['meter_reading'], label = 'Validation')\nplt.plot(pred_bidX_meterY['pred_meter_reading'], label = 'Holt Model')\nplt.legend(loc = 'best')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Number time!"},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Holt - RMSLE value:\")\nprint(rmsle(valid_holt[\"pred_meter_reading\"],\n           valid_holt[\"cur_meter_reading\"]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"holt_rmsle_list=[]\nfor meter_t in list(valid_holt[\"meter\"].unique()):\n        sub_valid_holt_df = valid_holt.loc[(\n            valid_holt[\"meter\"]==meter_t),:].copy()\n        sub_rmsle = rmsle(sub_valid_holt_df[\"pred_meter_reading\"],\n           sub_valid_holt_df[\"cur_meter_reading\"])\n        sub_rmsle_df = pd.DataFrame({\"meter\":[meter_t],\n                                   \"rmsle\":[sub_rmsle]})\n        holt_rmsle_list.append(sub_rmsle_df)\nholt_rmsle_list = pd.concat(holt_rmsle_list)\nholt_rmsle_list  ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Winter Holt Model\n **Winter Holt's model**  takes both trend and seasonality into account. Seasonality is found when time series data show variations a specific regular intervals less than a year (eg. hourly, weekly, monthly, quarterly). Unfortunately, Quarterly may be difficult to see with only a years worth of data, but keep an eyes out for other patterns in your EDA.  "},{"metadata":{},"cell_type":"markdown","source":"Here, since we split our training and validation set into 9 months in the training set and 3 months in the validation set. "},{"metadata":{"trusted":true},"cell_type":"code","source":"# Winter-Holt's prediction model\n# parameters\n# - train_dataframe: dataframe containing training data\n# - timeIdx_train: data frame with time series in the index\n# - valid_winterHolt: copy of the validation data frame\n# - seasonility: a list of all the seasonal period that are being tested\n# - b_id: building ID\n# - meter_t: meter number\n# - pred_col: list of the predicion column names (must be unique)\n# - plot: set to True to print a plot of the predictions\ndef winHolt(train_dataframe, timeIdx_train, valid_winterHolt, \n            seasonality, b_id, meter_t, plot=False, \n            pred_col=[], known_true=False):\n    if len(seasonality) > len(pred_col):\n        if len(seasonality) ==1:\n            pred_col=[\"winterHolt\"]\n        else:\n            pred_col=[]\n            for i in range(0,len(seasonality)):\n                pred_col.append(\"winterHolt_sp_\"+str(seasonality[i]))\n    ignored_pred_cols=[]\n    if(not ((b_id in train_dataframe[\"building_id\"]) and\n           (meter_t in train_dataframe.loc[train_dataframe[\"building_id\"]==b_id,\"meter\"].values))):\n        print(\"missing!\")\n        print(b_id)\n        print(meter_t)\n        # if there is no meter reading for a specific\n        # building ID then I'll just train\n        # independent of the building ID\n        sub_timeTrain_df = timeIdx_train.loc[(\n            timeIdx_train[\"meter\"]==meter_t),\"meter_reading\"].copy()\n        numValid = len(valid_winterHolt.loc[(\n            (valid_winterHolt[\"building_id\"]==b_id) &\n            (valid_winterHolt[\"meter\"]==meter_t)),:])\n        for i in range(0,len(seasonality)):\n            if (len(sub_timeTrain_df.index.unique())/\n                 seasonality[i]) >= 2:\n                fit_wintHolt = ExponentialSmoothing(\n                    sub_timeTrain_df,\n                    seasonal_periods=seasonality[i],\n                    trend='add',\n                    seasonal='add').fit()\n                # forecast the meter_readings\n                valid_winterHolt.loc[(\n                    (valid_winterHolt[\"building_id\"]==b_id) &\n                    (valid_winterHolt[\"meter\"]==meter_t)),pred_col[i]] = \\\n                    fit_wintHolt.forecast(numValid).values\n            else:\n                ignored_pred_cols.append(pred_col[i])\n    else:\n        # fit the model to the meter values of\n        # this building type\n        sub_timeTrain_df = timeIdx_train.loc[(\n                        (timeIdx_train[\"building_id\"]==b_id) &\n                        (timeIdx_train[\"meter\"]==meter_t)),\"meter_reading\"].copy()\n        numValid = len(valid_winterHolt.loc[(\n                        (valid_winterHolt[\"building_id\"]==b_id) &\n                        (valid_winterHolt[\"meter\"]==meter_t)),:])\n        for i in range(0,len(seasonality)):\n            if (len(sub_timeTrain_df.index.unique())/\n                 seasonality[i]) >= 2:\n                fit_wintHolt = ExponentialSmoothing(\n                    sub_timeTrain_df,\n                    seasonal_periods=seasonality[i],\n                    trend='add',\n                    seasonal='add').fit()\n                # forecast the meter_readings\n                valid_winterHolt.loc[(\n                    (valid_winterHolt[\"building_id\"]==b_id) &\n                    (valid_winterHolt[\"meter\"]==meter_t)),pred_col[i]] = \\\n                    fit_wintHolt.forecast(numValid).values\n            else:\n                ignored_pred_cols.append(pred_col[i])\n    if plot:\n        b_i=b_id\n        m_t=meter_t\n        train_bidX_meterY = train_dataframe.loc[(\n            (train_dataframe[\"building_id\"]==b_i) &\n            (train_dataframe[\"meter\"]==m_t)),:].copy()\n        \n        valid_bidX_meterY = valid_winterHolt.loc[(\n            (valid_winterHolt[\"building_id\"]==b_i) &\n            (valid_winterHolt[\"meter\"]==m_t)),:].copy()\n        plt.figure(figsize =(15,8))\n        plt.plot(train_bidX_meterY['meter_reading'], label = 'Train')\n        if known_true:\n            plt.plot(valid_bidX_meterY['cur_meter_reading'], label = 'Validation')\n        for i in range(0,len(seasonality)):\n            if pred_col[i] not in ignored_pred_cols:\n                plt.plot(valid_bidX_meterY[pred_col[i]], label = pred_col[i])\n        plt.legend(loc = 'best')\n    return(valid_winterHolt)    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# create new data frame for this model\nvalid_winterHolt = valid.copy()\n# rename timestamp to signify the current meter reading time\nvalid_winterHolt = valid_winterHolt.rename(\n    columns={\"timestamp\": \"now\", \n             \"meter_reading\": \"cur_meter_reading\",\n            \"log_meter_reading\":\"cur_log_meter_reading\"})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# set seasonality\nsp=[365-90,4*9,9,3]\n\n\n# This model splits the data based on \n# building ID and model type\nnbuildings=len(valid[\"building_id\"].unique())\nprint(\"number of buildings: \"+ str(nbuildings))\nx=0\nfor b_id in list(valid[\"building_id\"].unique()):\n    \n    update_progress(x / nbuildings)\n    x+=1\n    for meter_t in list(\n        valid_winterHolt.loc[valid_winterHolt[\"building_id\"]==b_id,\"meter\"].unique()):\n        print(b_id)\n        print(meter_t)\n        valid_winterHolt = winHolt(train,\n            timeIndexed_train, valid_winterHolt,\n            sp, b_id,meter_t)\nupdate_progress(1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's see how it looks compared to the true values in meters of some of the buildings"},{"metadata":{"trusted":true},"cell_type":"code","source":"sp=[365-90,4*9,9,3]\nb_id=0\nmeter_t=0\nvalid_winterHolt = winHolt(train, timeIndexed_train, valid_winterHolt, sp, b_id,meter_t, True, known_true=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sp=[365-90,4*9,9,3]\nb_id=161\nmeter_t=1\nvalid_winterHolt = winHolt(train, timeIndexed_train, valid_winterHolt, sp, b_id,meter_t, True, known_true=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sp=[365-90,4*9,9,3]\nb_id=745\nmeter_t=2\nvalid_winterHolt = winHolt(train, timeIndexed_train, valid_winterHolt, sp, b_id,meter_t, True, known_true=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sp=[365-90,4*9,9,3]\nb_id=106\nmeter_t=3\nvalid_winterHolt = winHolt(train, timeIndexed_train, valid_winterHolt, sp, b_id,meter_t, True, known_true=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Number time!"},{"metadata":{"trusted":true},"cell_type":"code","source":"pred_col=[]\nfor i in sp:\n    pred_colName = \"winterHolt_sp_\"+str(i)\n    print(\"winterHolt (sp =\"+str(i)+\") - RMSLE value:\")\n    print(rmsle(valid_winterHolt[pred_colName],\n           valid_winterHolt[\"cur_meter_reading\"]))\n    pred_col.append(pred_colName)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"winterHolt_rmsle_list=[]\nfor i in range(0,len(sp)):\n    for meter_t in list(valid_winterHolt[\"meter\"].unique()):\n        sub_valid_winterHolt_df = valid_winterHolt.loc[(\n            valid_winterHolt[\"meter\"]==meter_t),:].copy()\n        sub_rmsle = rmsle(sub_valid_winterHolt_df[pred_col[i]],\n           sub_valid_winterHolt_df[\"cur_meter_reading\"])\n        sub_rmsle_df = pd.DataFrame({\"seasonality\":[sp[i]],\n                                     \"meter\":[meter_t],\n                                   \"rmsle\":[sub_rmsle]})\n        winterHolt_rmsle_list.append(sub_rmsle_df)\nwinterHolt_rmsle_list = pd.concat(winterHolt_rmsle_list)\nprint(winterHolt_rmsle_list)  ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Summary"},{"metadata":{},"cell_type":"markdown","source":"To summarize our findings for each meter:\n* we predicted meter0 (electricity) with the smallest error rate (rmsle~1557) using Holt's model which ignores seasonality\n* we predicted meter1 (chilledwater) best using Winter-Holt's model and weekly smoothing (seasonality=36, rmsle~2045) or seasonal smoothing (seasonality=3, rmsle~2049)\n* meter2 (steam) works best by Winter Holt's model and daily smoothing (seasonality=275, rmsle~1390)\n* meter3 (hotwater) works best by Winter Holt's model and seasonal smoothing (seasonality=3, rmsle~1117)"},{"metadata":{},"cell_type":"markdown","source":"# Work I Learned From\n* https://www.youtube.com/watch?v=d4Sn6ny_5LI (Time Series Forecasting Video)"}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":1}