{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Time series modeling\n", "\n", "files needed = ( none )\n", "\n", "We continue to work with the statsmodels package [(docs)](https://devdocs.io/statsmodels/). This notebook is a brief introduction to the time series methods. \n", "\n", "We are only going to touch on autoregressive models here, but statsmodels [tsa](https://www.statsmodels.org/dev/tsa.html) package has support for moving average and ARIMA models as well as vector autoregression and error correction models. The package also includes filtering methods and time series related plotting functions. \n", "\n", "Time series models are often used for forecasting. If forecasting interests you, Prof. Hansen offers a [forecasting course](https://www.ssc.wisc.edu/~bhansen/460/) (ECON 460) where you can learn a lot more about these kinds of models." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import pandas as pd # for data handling\n", "import numpy as np # for numerical methods and data structures\n", "import matplotlib.pyplot as plt # for plotting\n", "import seaborn as sea # advanced plotting\n", "\n", "import statsmodels.graphics.tsaplots as tsaplots # Gives the the autocorrelation plot\n", "import statsmodels.tsa as tsa # The time series models\n", "\n", "from pandas_datareader import data, wb # we are grabbing the data and wb functions from the package\n", "import datetime as dt # for time and date" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### GDP \n", "\n", "Let's model US GDP as an autoregressive process. Grab the data from FRED." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "codes = ['GDPC1'] # The code for real GDP at FRED.\n", "start = start = dt.datetime(1947, 1, 1)\n", "\n", "usa = data.DataReader(codes, 'fred', start)\n", "usa.info()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "usa.rename(columns={'GDPC1':'gdp'}, inplace=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(15,6))\n", "\n", "ax.plot(usa.index, usa['gdp'], color = 'red')\n", "ax.set_ylabel('bil 2012 dollars')\n", "\n", "sea.despine(ax=ax)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Stationarity\n", "\n", "The GDP data are not stationary --- there is clearly a trend. Autoregressive models are usually used on stationary data. We can filter the data to remove the trend and recover a stationary series (tsa has [methods](https://www.statsmodels.org/dev/tsa.html#time-series-filters) to to this), but GDP is usually stationary in growth rates. \n", "\n", "Compute the growth rate of the GDP data." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# The log-difference is a close approximation of the growth rate when growth rates are small. \n", "# log-differences are often used in practice because they are symmetric. \n", "\n", "usa['gdp_diff'] = np.log(usa['gdp']) - np.log(usa['gdp'].shift(1))\n", "usa['gdp_pct'] = usa['gdp'].pct_change()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(15,6))\n", "\n", "ax.plot(usa.index, usa['gdp_diff'], color = 'red')\n", "ax.plot(usa.index, usa['gdp_pct'], color = 'blue')\n", "\n", "ax.set_ylabel('growth rate')\n", "\n", "sea.despine(ax=ax)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This looks roughly stationary. There are things like Dickey-Fuller tests for stationarity (again, these can be done in statsmodels) but they are outside of the scope of this notebook. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Autoregressive model \n", "\n", "We say [autoregressive model](https://en.wikipedia.org/wiki/Autoregressive_model) because we are modeling the variable as a function of its past values. The general model is\n", "\n", "$$y_t = \\varphi_0 + \\varphi_1y_{t-1} + \\varphi_2 y_{t-2} + \\cdots + \\varphi_p y_{t-p}+\\epsilon_t.$$\n", "\n", "This model has $p$ lags. We refer to it as an AR(p) model or an autoregressive model of order $p$. \n", "\n", "Determining the number of lags to include is part of specifying the model.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's start by looking at the data to see if it has an autoregressive property. Is there a relationship between the current value of GDP growth and its past value?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(10,10))\n", "\n", "ax.scatter(usa['gdp_diff'], usa['gdp_diff'].shift(1), color = 'red')\n", "\n", "# plot the 45 degree line\n", "ax.plot([-0.02, 0.045], [-0.02, 0.045], color='black')\n", "\n", "ax.text(0.04, 0.037, '45-degree line')\n", "\n", "ax.set_ylabel('gdp(t)')\n", "ax.set_xlabel('gdp(t-1)')\n", "\n", "sea.despine(ax=ax, trim=True)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There appears to be a positive relationship, but it is noisy. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The autocorrelation function\n", "\n", "Our plot above tells us about the $t$ and $t-1$ relationship. How about $t$ and $t-2$? $t-5$? We could continue to make the plots above, but there are better ways. \n", "\n", "The autocorrelation function plots $\\text{corr}(y_t, y_{t-k})$ for many values of $k$. statsmodels automates this for us with the .plot_acf() function [(docs)](https://www.statsmodels.org/dev/generated/statsmodels.graphics.tsaplots.plot_acf.html). " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(10,6))\n", "\n", "# plot_acf is picky about nas, so drop them\n", "# The 'lags' parameter determines how many lags to show\n", "tsaplots.plot_acf(usa['gdp_diff'].dropna(), lags = 8, ax = ax)\n", "\n", "ax.set_xlabel('lag = k', fontsize=14)\n", "ax.set_ylabel(r'corr(gdp$_t$, gdp$_{t-k})$', fontsize=14)\n", "sea.despine(ax=ax)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Estimating the AR(p) model\n", "\n", "The AR(p) model in statsmodels works like the other models. You create a model object, then you call the .fit( ) method to estimate the parameters. \n", "\n", "The only real difference here is that you do not need to write out a string with the model specification. The .AR( ) method knows that the model is AR(p). " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# I couldn't get these two methods to work. Need to look into this.\n", "# usa = usa.asfreq('q') \n", "# usa.index.freq = 'q'\n", "\n", "# Need to set the frequency of hte data. This is a not-very-elegant way. \n", "usa = usa.resample('q').mean()\n", "\n", "# Construct the model. \n", "ar_mod = tsa.ar_model.AR(usa['gdp_diff'].dropna(), freq='Q')\n", "\n", "type(ar_mod) # What have we got?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we call fit [(docs)](https://www.statsmodels.org/dev/generated/statsmodels.tsa.ar_model.AR.fit.html#statsmodels.tsa.ar_model.AR.fit). \n", "\n", "Here we have to decide how many lags we want. There are two options to deal with. \n", "\n", "* ic is the criterion for choosing how many lags to include. Options are 'aic', 'bic', 'hic', and 't-stat'.Which one to choose is a technical detail. Let's go with 'aic'. Alternatively, you can choose ic = None.\n", "* maxlag is the maximum number of lags. If ic = None, this is the number of lags. \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#res = ar_mod.fit(ic = None, maxlag=3)\n", "\n", "res = ar_mod.fit(ic = 'aic')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# The parameters\n", "print(res.params)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# The optimal number of lags is\n", "print('The number of lags is {0}.'.format(res.k_ar))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The results object has lots of stuff in it. Try tab completion to see what's in there. \n", "\n", "Let's plot the fitted values" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(10,6))\n", "\n", "ax.plot(res.fittedvalues, color='blue', label = 'fitted values')\n", "ax.plot(usa.index, usa.gdp_diff, color='black', alpha = 0.5, label = 'data')\n", "\n", "sea.despine(ax=ax)\n", "ax.set_ylabel('growth rate')\n", "ax.set_title('US GDP growth rates')\n", "\n", "ax.legend(frameon=False)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The fitted values track the data, but do not generate enough volatility. If we really wanted to get serious about modeling GDP growth, we would want to expand our model to something like an autoregressive-moving-average model." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Practice\n", "\n", "Take a few minutes and try the following. Feel free to chat with those around if you get stuck. The TA and I are here, too." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. Read in the S\\&P 500 price index from FRED. Use November 23, 2013 as the start date.\n", "2. Resample the data to weekly frequency using mean()." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3. Plot the weekly average price. Are the data stationary?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "4. Compute the log difference of the weekly price data. \n", "5. Plot the growth rates. Are the data stationary?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "6. Plot the autocorrelation function for 20 lags. What lags look meaningful? " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "7. Estimate an AR(p) model of the price data. Include 1 lag." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "8. Print out the number of lags used and the parameter estimates. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "9. Plot the data and the fitted values. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "10. Restimate the model letting the .fit( ) method choose the number of lags. How many lags were chosen?\n", "11. Print out the number of lags and the parameter estimates. \n", "\n", "\$If you have extra time, try ploting the estimates and their confidence intervals. I used errorbar, but there are many ways.\$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "12. Plot the data and the fitted values from the 1-lag model and the 2-lag model. Plot just for the period 1/1/2015 through 1/1/2017." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Attachments", "kernelspec": { "display_name": "Python 3", "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.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }