{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install causal-impact","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Python Libraries to know\n- Pandas: For loading and manipulating data. \n- NumPy: Working with numbers and matrices. \n- Matplotlib (Seaborn too): For plotting and graphing data. Matplotlib will do the basics, and Seaborn can do some more advanced plots. \n- Sklearn: Has a bunch of different models and useful functions in an easy to use format.\n- Statsmodels: Many models, time series models and statistical tests. \n\n### Nice to know (but not necessary)\n- PyTorch or Tensorflow: Neural networks and gradient descent optimization. Last year there was text data for both competitions I went to, but I'm under the impression the winners didn't rely heavily on NLP.\n- LightGBM or XGBoost: Easy to use black-box gradient boosting trees models. Handles missing values automatically. You can use shapley values to \"interepret\" the model. \n- PyMC3, PyStan or other PPL: Very powerful but has a learning curve. For bayesian modeling.  \n- CVXPY: Convex optimization. Useful if you are trying to optimize over a convex function. You can always use a general optimizer (ie particle swarm, bayesian optimization, grid-search) instead. \n\n### Useful Resources:\nKaggle\n\nhttps://datasetsearch.research.google.com/"},{"metadata":{"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\ntrain_df = pd.read_csv('../input/ashrae-energy-prediction/train.csv')\nweather_df = pd.read_csv('../input/ashrae-energy-prediction/weather_train.csv')\nmeta_df = pd.read_csv('../input/ashrae-energy-prediction/building_metadata.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"##  Look at what kind of data is in these files"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_df.head()\n# There are some missing values. We should also eventually ensure that all of the values fall within a reasonable range. ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"meta_df.head()\n# Missing values as well. ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Explore train_df"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df['meter'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df['timestamp'][0] ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Timestamp is not in a date time format, so let's convert to the pandas date format, and then add some additional features!"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df['timestamp'] = pd.to_datetime(train_df['timestamp'])\n\ntrain_df['month'] = train_df['timestamp'].dt.month\ntrain_df['weekday'] = train_df['timestamp'].dt.dayofweek\ntrain_df['monthday'] = train_df['timestamp'].dt.day\ntrain_df['hour'] = train_df['timestamp'].dt.hour\ntrain_df['minute'] = train_df['timestamp'].dt.minute","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df['minute'].unique() # Looks like the data doesn't go down to minute resolution. Lets drop it. ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = train_df.drop(['minute'], axis = 1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Look at individual buildings."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(train_df[train_df['building_id'] == 0]['meter_reading'], alpha = 0.8)\nplt.plot(train_df[train_df['building_id'] == 1]['meter_reading'], alpha = 0.8)\nplt.plot(train_df[train_df['building_id'] == 2]['meter_reading'], alpha = 0.8)\nplt.plot(train_df[train_df['building_id'] == 500]['meter_reading'], alpha = 0.8)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Let's look at the autocorrelation of these plots, and look at a lagplot."},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.plotting.lag_plot(train_df[train_df['building_id'] == 0]['meter_reading'])\nplt.plot([0,400],[0,400])\n# Look at the 3 clusters. ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.plotting.lag_plot(train_df[train_df['building_id'] == 500]['meter_reading'])\nplt.plot([0,400],[0,400])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.plotting.autocorrelation_plot(train_df[train_df['building_id'] == 500]['meter_reading'])\nplt.show()\npd.plotting.autocorrelation_plot(train_df[train_df['building_id'] == 500]['meter_reading'][:300])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Some buildings have multiple power meters."},{"metadata":{"scrolled":true,"trusted":true},"cell_type":"code","source":"train_df[train_df['meter'] == 2].head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### The meters are not necessarily consecutive numbers."},{"metadata":{"scrolled":true,"trusted":true},"cell_type":"code","source":"train_df[train_df['building_id'] == 745].head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(train_df[train_df['building_id'] == 745].meter.unique())\nprint(train_df[train_df['building_id'] == 1414].meter.unique())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### What does the distribution of weather look like?"},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(weather_df['air_temperature'].dropna())\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### We have 3 dataframes, but they should be merged into one so we can feed it to a model.\n### The site_id will be mapped to the building_id between the train and the meta, and the weather will be mapped to the site and time of the training data."},{"metadata":{"trusted":true},"cell_type":"code","source":"all_df = pd.merge(train_df, meta_df, on = 'building_id', how = 'left')\nall_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"weather_df['timestamp'] = pd.to_datetime(weather_df['timestamp']) # Convert weather to the correct format before merging\nall_df = pd.merge(all_df, weather_df, on = ['site_id', 'timestamp'], how = 'left')\nall_df['date'] = all_df['timestamp'].dt.date\nall_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Now we have the data in a format we can use, so we now need to think of some problems we can solve with this data. It is a good idea to Google what relevant news stories there are around power consumption, buildings, and weather. Then also look to see if there is any relevant research and papers about the topic. \n\n### Some ideas I came up with are:\n### 1. New York recently came up with a tax on inefficient buildings. How will this tax affect power consumption?\n### 2. How will climate change affect power consumption?\n### 3. What is the most effective way to transition to renewable power? What combinaiton of Solar, Wind and Batteries would be the most cost effective and would be robust to prolonged bad weather. \n"},{"metadata":{},"cell_type":"markdown","source":"\n### Then consider some relevant factors such as: How difficult it will be to answer the question, what external data can be brough in, is there any relevant research, how much will the topic impress judges?\n\n### 1. We would likely need some electricity/building pricing information or data before and after the tax/policy change to be able to solve it. This is an active topic of discussion in NYC, so there are probably many different opinions that can be discussed in good narrative. There might be difficulties since data is noisy, and there might not be a visible reaction from changes to electricity prices (risky). Maybe there are some other policies which can be better examined, but if this can be pulled off, it will be a top contender. \n### 2. This is very do-able, even with just linear regression. Not much external data required, just some estimates of the effects of climate change. Just create a model which predicts power usage from weather and adjust the weather to climate change predictions. Might be less impressive to the judges, but has serious potential to win if done right.\n### 3.  Solvable, but harder than 2. Would need pricing data on renewable power, and how to convert weather to power generation. Very impressive if you manage to solve it (similar idea won us the Championship).\n\n### Let's take a crack at #2 since it is the easiest. "},{"metadata":{},"cell_type":"markdown","source":"### Causation\nCounterfactuals are one of the easiest ways to show causation. The goal of a counterfactual approach is to estimate how **Y** (energy use) would be different had **X** (climate) been something else. There are two common ways of establishing a counterfactual with data. \n\n- Matching: \nMatching based approaches try to replicate a controlled experiment. We can match buildings in cold weather to similar buildings in warm weather and compare differences in energy use between similar buildings. \nhttps://github.com/benmiroglio/pymatch\n- Model Based Counterfactuals:\nThe idea here is to \"learn\" a counter factual with a model. One way we can do this is to create a forecast assuming the temperature stays the same, and compare that forecast to the real results once the temperature changes. This isn't a great example since the weather is always changing and this approach works best with point-in-time interventions, but it illustrates how to use this approach.\nhttps://github.com/tcassou/causal_impact"},{"metadata":{"trusted":true},"cell_type":"code","source":"data = all_df.groupby(['date', 'monthday'])[['meter']].mean() \\\n        .join(all_df.groupby(['date', 'monthday'])['air_temperature'].mean()).sort_values(['date']).reset_index()\ndata.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data['air_temperature'].iloc[200:500].plot()\nplt.plot([264,264],[14,28])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from causal_impact import CausalImpact\ndata['weekend'] = (pd.to_datetime(data.date).dt.dayofweek >= 5).astype(int) / 100 + 0.66\ndata_time = data[['meter','weekend']].iloc[200:350].rename({'meter':'y','weekend':'x1'}, axis = 1).reset_index(drop = True)\n\nci = CausalImpact(data_time, 264-200)\nci.run(max_iter=1000)\nci.plot()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data['meter'].iloc[200:500].plot()\nplt.plot([264,264],[0.66,0.68])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from statsmodels.regression.linear_model import OLS\ndata['constant'] = 1\nmodel = OLS(data['meter'], data[['air_temperature','weekend', 'constant']])\nresults = model.fit()\nresults.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(results.resid)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(results.resid)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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":4}