{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from random import gauss\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from arch import arch_model\n",
    "from statsmodels.graphics.tsaplots import plot_acf, plot_pacf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GARCH(2,2) Model\n",
    "\n",
    "$$\n",
    "a_t = \\varepsilon_t \\sqrt{\\omega + \\alpha_1 a_{t-1}^2 + \\alpha_2 a_{t-2}^2 + \\beta_1 \\sigma_{t-1}^2 + \\beta_2 \\sigma_{t-2}^2}\n",
    "$$\n",
    "\n",
    "$$\n",
    "a_0, a_1 \\sim \\mathcal{N}(0,1)\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\sigma_0 =1, \\sigma_1 = 1\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\varepsilon_t \\sim \\mathcal{N}(0,1)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create dataset\n",
    "n = 1000\n",
    "omega = 0.5\n",
    "\n",
    "alpha_1 = 0.1\n",
    "alpha_2 = 0.2\n",
    "\n",
    "beta_1 = 0.3\n",
    "beta_2 = 0.4\n",
    "\n",
    "test_size = int(n*0.1)\n",
    "\n",
    "series = [gauss(0,1), gauss(0,1)]\n",
    "vols = [1, 1]\n",
    "\n",
    "for _ in range(n):\n",
    "    new_vol = np.sqrt(omega + alpha_1*series[-1]**2 + alpha_2*series[-2]**2 + beta_1*vols[-1]**2 + beta_2*vols[-2]**2)\n",
    "    new_val = gauss(0,1) * new_vol\n",
    "    \n",
    "    vols.append(new_vol)\n",
    "    series.append(new_val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,4))\n",
    "plt.plot(series)\n",
    "plt.title('Simulated GARCH(2,2) Data', fontsize=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,4))\n",
    "plt.plot(vols)\n",
    "plt.title('Data Volatility', fontsize=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,4))\n",
    "plt.plot(series)\n",
    "plt.plot(vols, color='red')\n",
    "plt.title('Data and Volatility', fontsize=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PACF Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_pacf(np.array(series)**2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fit the GARCH Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train, test = series[:-test_size], series[-test_size:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = arch_model(train, p=2, q=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_fit = model.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_fit.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Predict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "predictions = model_fit.forecast(horizon=test_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,4))\n",
    "true, = plt.plot(vols[-test_size:])\n",
    "preds, = plt.plot(np.sqrt(predictions.variance.values[-1, :]))\n",
    "plt.title('Volatility Prediction', fontsize=20)\n",
    "plt.legend(['True Volatility', 'Predicted Volatility'], fontsize=16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "predictions_long_term = model_fit.forecast(horizon=1000)\n",
    "plt.figure(figsize=(10,4))\n",
    "true, = plt.plot(vols[-test_size:])\n",
    "preds, = plt.plot(np.sqrt(predictions_long_term.variance.values[-1, :]))\n",
    "plt.title('Long Term Volatility Prediction', fontsize=20)\n",
    "plt.legend(['True Volatility', 'Predicted Volatility'], fontsize=16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rolling Forecast Origin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rolling_predictions = []\n",
    "for i in range(test_size):\n",
    "    train = series[:-(test_size-i)]\n",
    "    model = arch_model(train, p=2, q=2)\n",
    "    model_fit = model.fit(disp='off')\n",
    "    pred = model_fit.forecast(horizon=1)\n",
    "    rolling_predictions.append(np.sqrt(pred.variance.values[-1,:][0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,4))\n",
    "true, = plt.plot(vols[-test_size:])\n",
    "preds, = plt.plot(rolling_predictions)\n",
    "plt.title('Volatility Prediction - Rolling Forecast', fontsize=20)\n",
    "plt.legend(['True Volatility', 'Predicted Volatility'], fontsize=16)"
   ]
  }
 ],
 "metadata": {
  "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
