arch_testing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import arch\n",
    "from arch.univariate import Normal as ARCHNormal\n",
    "# np.set_printoptions(precision=2, suppress=True, linewidth=100, floatmode='fixed')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set seed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = 42"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Random number generator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\sigma_t^2 = \\omega + \\alpha \\epsilon_{t-1}^2$  \n",
    "Below we have $y_t = 1$ for all $t$, a zero mean model and $\\omega = 0$ and $\\alpha = 1$ so $\\sigma_t^2 = 1$ for all $t$  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.305, -1.04 ,  0.75 ,  0.941, -1.951, -1.302,  0.128, -0.316,\n",
       "       -0.017, -0.853])"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# get the generated normal random numbers from the arch.arch_model.forecast function with the same seed\n",
    "length = 10\n",
    "test_df = pd.DataFrame(np.ones(length), columns=['y'])\n",
    "test_model = arch.arch_model(y=test_df['y'], x=None, mean='Zero', vol='ARCH', p=1, rescale=False)\n",
    "test_model.distribution = ARCHNormal(seed=seed)\n",
    "sims = test_model.forecast(params=[0,1], horizon=1, start=0, method='simulation', simulations=1, x=None, align='origin')\n",
    "normal_rvs = sims.simulations.values.squeeze()\n",
    "normal_rvs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test vol='ARCH'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# create artificial data to test the arch.arch_model.forecast function\n",
    "# data will contain integer values so that it is easier to compare the results\n",
    "test_df = pd.DataFrame(np.arange(10), columns=['y'])\n",
    "y = test_df['y'].values\n",
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration:      1,   Func. Count:      4,   Neg. LLF: 28.43922788393828\n",
      "Iteration:      2,   Func. Count:      7,   Neg. LLF: 28.427130663050598\n",
      "Iteration:      3,   Func. Count:     10,   Neg. LLF: 28.423630703038214\n",
      "Iteration:      4,   Func. Count:     13,   Neg. LLF: 28.42333130448614\n",
      "Iteration:      5,   Func. Count:     16,   Neg. LLF: 28.423323726474315\n",
      "Iteration:      6,   Func. Count:     18,   Neg. LLF: 28.42332373612137\n",
      "Optimization terminated successfully    (Exit mode 0)\n",
      "            Current function value: 28.423323726474315\n",
      "            Iterations: 6\n",
      "            Function evaluations: 18\n",
      "            Gradient evaluations: 6\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "                        Zero Mean - ARCH Model Results                        \n",
       "==============================================================================\n",
       "Dep. Variable:                      y   R-squared:                       0.000\n",
       "Mean Model:                 Zero Mean   Adj. R-squared:                  0.100\n",
       "Vol Model:                       ARCH   Log-Likelihood:               -28.4233\n",
       "Distribution:                  Normal   AIC:                           60.8466\n",
       "Method:            Maximum Likelihood   BIC:                           61.4518\n",
       "                                        No. Observations:                   10\n",
       "Date:                Thu, Oct 31 2024   Df Residuals:                       10\n",
       "Time:                        07:10:30   Df Model:                            0\n",
       "                            Volatility Model                            \n",
       "========================================================================\n",
       "                 coef    std err          t      P>|t|  95.0% Conf. Int.\n",
       "------------------------------------------------------------------------\n",
       "omega          1.9896      1.452      1.370      0.171 [ -0.856,  4.835]\n",
       "alpha[1]       1.0000      0.196      5.106  3.283e-07 [  0.616,  1.384]\n",
       "========================================================================\n",
       "\n",
       "Covariance estimator: robust\n",
       "ARCHModelResult, id: 0x12f3ef3d0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# only one lag on residual with vol='ARCH'\n",
    "test_model = arch.arch_model(y=test_df['y'], x=None, mean='Zero', vol='ARCH', p=1, rescale=False)\n",
    "test_model.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1\n",
    "$\\sigma_t^2 = \\omega + \\alpha_1 \\epsilon_{t-1}^2$  \n",
    "Below we have $y_t = t$ for $t=0,\\ldots,9$, a zero mean model, $\\omega = 0$ and $\\alpha_1 = 1$  \n",
    "So $\\sigma_t^2 = y_{t-1}^2$ for $t=1,\\ldots,10$  \n",
    "e.g. $\\sigma_5^2 = y_4^2 = 4^2 = 16$  \n",
    "Note that vars[i] $= \\text{VAR}[\\hat{y}_{i+1}]$ for $i=0,\\ldots,9$ since start=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0.  1.  4.  9. 16. 25. 36. 49. 64. 81.]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_model.distribution = ARCHNormal(seed=seed)\n",
    "sims = test_model.forecast(params=[0,1], horizon=1, start=0, method='simulation', simulations=1, x=None, align='origin')\n",
    "series = sims.simulations.values.squeeze()\n",
    "means = sims.mean.values.squeeze()\n",
    "vars = sims.variance.values.squeeze()\n",
    "print(means, vars)\n",
    "np.allclose(vars, y**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2\n",
    "$\\sigma_t^2 = \\omega + \\alpha_1 \\epsilon_{t-1}^2 + \\alpha_2 \\epsilon_{t-2}^2$  \n",
    "Below we have $y_t = t$ for $t=0,\\ldots,9$, a zero mean model, $\\omega = 0$ and $\\alpha_1 = \\alpha_2 = 1$  \n",
    "So $\\sigma_t^2 = y_{t-1}^2 + y_{t-2}^2$ for $t=1,\\ldots,10$  \n",
    "e.g. $\\sigma_5^2 = y_4^2 + y_3^2 = 4^2 + 3^2 = 16 + 9 = 25$  \n",
    "Note that vars[i] $= \\text{VAR}[\\hat{y}_{i+2}]$ for $i=0,\\ldots,8$ since start=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 0. 0. 0. 0. 0. 0. 0. 0.] [  1.   5.  13.  25.  41.  61.  85. 113. 145.]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# two lags on residual with vol='ARCH' and start=1\n",
    "test_model = arch.arch_model(y=test_df['y'], x=None, mean='Zero', vol='ARCH', lags=0, p=2, rescale=False)\n",
    "test_model.distribution = ARCHNormal(seed=seed)\n",
    "sims = test_model.forecast(params=[0,1,1], horizon=1, start=1, method='simulation', simulations=2, x=None, align='origin')\n",
    "series = sims.simulations.values.squeeze()\n",
    "# print(f'series.shape: {series.shape}')\n",
    "means = sims.mean.values.squeeze()\n",
    "vars = sims.variance.values.squeeze()\n",
    "print(means, vars)\n",
    "np.allclose(vars, y[1:]**2 + y[:-1]**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 3\n",
    "When start=None (default), the forecast starts from the last available observation.  \n",
    "Therefore, vars[0] $= \\text{VAR}[\\hat{y}_{10}] = 9^2 + 8^2 = 145$ as in Example 2  \n",
    "vars[1] $ = \\text{VAR}[\\hat{y}_{11}] = y_{10}^2 + y_{9}^2$  \n",
    "vars[2] $ = \\text{VAR}[\\hat{y}_{12}] = y_{11}^2 + y_{10}^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "values: [  3.66927954 -10.10785438   8.06978747]\n",
      "means: [0. 0. 0.]\n",
      "vars: [145.          94.46361231 115.63233254]\n",
      "y_9^2 + y_8^2: 94.46361231063993\n",
      "y_10^2 + y_9^2: 115.63233253526381\n"
     ]
    }
   ],
   "source": [
    "# two lags on residual with vol='ARCH' and start=None\n",
    "test_model = arch.arch_model(y=test_df['y'], x=None, mean='Zero', vol='ARCH', lags=0, p=2, rescale=False)\n",
    "test_model.distribution = ARCHNormal(seed=seed)\n",
    "sims = test_model.forecast(params=[0,1,1], horizon=3, method='simulation', simulations=1, x=None, align='origin')\n",
    "values = sims.simulations.values.squeeze()\n",
    "means = sims.mean.values.squeeze()\n",
    "vars = sims.variance.values.squeeze()\n",
    "print(f'values: {values}')\n",
    "print(f'means: {means}')\n",
    "print(f'vars: {vars}')\n",
    "print(f'y_9^2 + y_8^2: {values[0]**2 + 9**2}')\n",
    "print(f'y_10^2 + y_9^2: {values[1]**2 + values[0]**2}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Show vol='ARCH' == vol='GARCH' with q=0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 4\n",
    "Below is the same as Example 1 but with vol='GARCH' and q=0  \n",
    "Therefore values are expected to be the same as in Example 1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.  1.  4.  9. 16. 25. 36. 49. 64. 81.]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# one lag on residual with vol='GARCH'\n",
    "test_model = arch.arch_model(y=test_df['y'], x=None, mean='Zero', vol='GARCH', lags=1, p=1, q=0, rescale=False)\n",
    "test_model.distribution = ARCHNormal(seed=seed)\n",
    "sims = test_model.forecast(params=[0,1], horizon=1, start=0, method='simulation', simulations=1, x=None, align='origin')\n",
    "series = sims.simulations.values.squeeze()\n",
    "vars = sims.variance.values.squeeze()\n",
    "print(vars)\n",
    "np.allclose(vars, y**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 5\n",
    "Below is the same as Example 2 but with vol='GARCH' and q=0  \n",
    "Therefore values are expected to be the same as in Example 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  1.   5.  13.  25.  41.  61.  85. 113. 145.]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# two lags on residual with vol='GARCH' and start=1\n",
    "test_model = arch.arch_model(y=test_df['y'], x=None, mean='Zero', vol='GARCH', lags=1, p=2, q=0, rescale=False)\n",
    "test_model.distribution = ARCHNormal(seed=seed)\n",
    "sims = test_model.forecast(params=[0,1,1], horizon=1, start=1, method='simulation', simulations=1, x=None, align='origin')\n",
    "series = sims.simulations.values.squeeze()\n",
    "vars = sims.variance.values.squeeze()\n",
    "print(vars)\n",
    "np.allclose(vars, y[1:]**2 + y[:-1]**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Start=0 with insufficient lags"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 6\n",
    "\n",
    "**<span style='color:red'>When there are insufficient lags to compute the variance, the code is drawing lag residuals from unknown source hence the values should not be used.</span>**  \n",
    "For unknown reasons, the lag residuals that are not supplied by the user appear to be the same number for all unavailable lags as shown below.   \n",
    "Constant value of (24.03524054)**0.5 is used regardless of the seed. It is unclear where this value is coming from."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[24.03524054  0.          1.          4.          9.         16.\n",
      " 25.         36.         49.         64.        ]\n"
     ]
    }
   ],
   "source": [
    "# coefficient is one only for t-2 residual\n",
    "test_model = arch.arch_model(y=test_df['y'], x=None, mean='Zero', vol='ARCH', lags=1, p=4, rescale=False)\n",
    "test_model.distribution = ARCHNormal(seed=seed)\n",
    "sims = test_model.forecast(params=[0,0,1,0,0], horizon=1, start=0, method='simulation', simulations=1, x=None, align='origin')\n",
    "series = sims.simulations.values.squeeze()\n",
    "vars = sims.variance.values.squeeze()\n",
    "print(vars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[24.03524054 24.03524054  0.          1.          4.          9.\n",
      " 16.         25.         36.         49.        ]\n"
     ]
    }
   ],
   "source": [
    "# coefficient is one only for t-3 residual\n",
    "test_model = arch.arch_model(y=test_df['y'], x=None, mean='Zero', vol='ARCH', lags=1, p=4, rescale=False)\n",
    "test_model.distribution = ARCHNormal(seed=seed)\n",
    "sims = test_model.forecast(params=[0,0,0,1,0], horizon=1, start=0, method='simulation', simulations=1, x=None, align='origin')\n",
    "series = sims.simulations.values.squeeze()\n",
    "vars = sims.variance.values.squeeze()\n",
    "print(vars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[24.03524054 24.03524054 24.03524054  0.          1.          4.\n",
      "  9.         16.         25.         36.        ]\n"
     ]
    }
   ],
   "source": [
    "# coefficient is one only for t-4 residual\n",
    "test_model = arch.arch_model(y=test_df['y'], x=None, mean='Zero', vol='ARCH', lags=1, p=4, rescale=False)\n",
    "test_model.distribution = ARCHNormal(seed=seed)\n",
    "sims = test_model.forecast(params=[0,0,0,0,1], horizon=1, start=0, method='simulation', simulations=1, x=None, align='origin')\n",
    "series = sims.simulations.values.squeeze()\n",
    "vars = sims.variance.values.squeeze()\n",
    "print(vars)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ARX model for mean with GARCH for volatility"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>y</th>\n",
       "      <th>x</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2.0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3.0</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>5.0</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>8.0</td>\n",
       "      <td>6.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>13.0</td>\n",
       "      <td>7.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>21.0</td>\n",
       "      <td>8.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>34.0</td>\n",
       "      <td>9.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>55.0</td>\n",
       "      <td>10.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       y     x\n",
       "0    0.0   0.0\n",
       "1    1.0   1.0\n",
       "2    1.0   2.0\n",
       "3    2.0   3.0\n",
       "4    3.0   4.0\n",
       "5    5.0   5.0\n",
       "6    8.0   6.0\n",
       "7   13.0   7.0\n",
       "8   21.0   8.0\n",
       "9   34.0   9.0\n",
       "10  55.0  10.0"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# create artificial data to test the arch.arch_model.forecast function\n",
    "# data will contain integer values so that it is easier to compare the results\n",
    "# test_df = pd.DataFrame({'y':np.arange(10)*0.2, 'x':np.arange(10)})\n",
    "test_df = pd.DataFrame({'y':np.array([0,1,1,2,3,5,8,13,21,34,55]).astype(float), 'x':np.arange(11).astype(float)})\n",
    "y = test_df['y'].values\n",
    "x = test_df.x.values\n",
    "test_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration:      1,   Func. Count:     23,   Neg. LLF: -2656.8083929322147\n",
      "Inequality constraints incompatible    (Exit mode 4)\n",
      "            Current function value: -2656.808392273039\n",
      "            Iterations: 1\n",
      "            Function evaluations: 23\n",
      "            Gradient evaluations: 1\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/luchungi/opt/anaconda3/envs/signature/lib/python3.11/site-packages/arch/univariate/base.py:766: ConvergenceWarning: The optimizer returned code 4. The message is:\n",
      "Inequality constraints incompatible\n",
      "See scipy.optimize.fmin_slsqp for code meaning.\n",
      "\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Const        0.0\n",
       "y[1]         1.0\n",
       "y[2]         0.0\n",
       "y[3]         0.0\n",
       "y[4]         0.0\n",
       "y[5]         0.0\n",
       "y[6]         0.0\n",
       "y[7]         0.0\n",
       "y[8]         0.0\n",
       "y[9]         0.0\n",
       "y[10]        0.0\n",
       "omega        0.0\n",
       "alpha[1]     1.0\n",
       "alpha[2]     0.0\n",
       "alpha[3]     0.0\n",
       "alpha[4]     0.0\n",
       "alpha[5]     0.0\n",
       "alpha[6]     0.0\n",
       "alpha[7]     0.0\n",
       "alpha[8]     0.0\n",
       "alpha[9]     0.0\n",
       "alpha[10]    0.0\n",
       "Name: params, dtype: float64"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "large_df = pd.DataFrame({'y':np.arange(100).astype(float)})\n",
    "test_model = arch.arch_model(y=large_df['y'], x=None, mean='ARX', vol='GARCH', lags=10, p=10, q=0, rescale=False)\n",
    "param_res = test_model.fit()\n",
    "param_res.params[:] = 0.\n",
    "param_res.params[1] = 1.\n",
    "param_res.params[12] = 1.\n",
    "param_res.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Need to pad the df with values else the test_model.fix() will fail due to insufficient data.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>y</th>\n",
       "      <th>x</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>106</th>\n",
       "      <td>8.0</td>\n",
       "      <td>6.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>107</th>\n",
       "      <td>13.0</td>\n",
       "      <td>7.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>108</th>\n",
       "      <td>21.0</td>\n",
       "      <td>8.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>109</th>\n",
       "      <td>34.0</td>\n",
       "      <td>9.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>110</th>\n",
       "      <td>55.0</td>\n",
       "      <td>10.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>111 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "        y     x\n",
       "0     0.0   NaN\n",
       "1     1.0   NaN\n",
       "2     2.0   NaN\n",
       "3     3.0   NaN\n",
       "4     4.0   NaN\n",
       "..    ...   ...\n",
       "106   8.0   6.0\n",
       "107  13.0   7.0\n",
       "108  21.0   8.0\n",
       "109  34.0   9.0\n",
       "110  55.0  10.0\n",
       "\n",
       "[111 rows x 2 columns]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "temp_df = pd.concat([large_df, test_df], ignore_index=True)\n",
    "temp_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let time $t$ be the last available observation.  \n",
    "$E[y_{t+1}] = y_t = 55.$  \n",
    "$\\text{VAR}[y_{t+1}] = \\epsilon_{t}^2 = (y_{t} - E[y_{t}])^2 = (55 - 34)^2 = 441$  \n",
    "$E[y_{t+2}] = y_{t+1} = $ values[0]  \n",
    "$\\text{VAR}[y_{t+2}] = \\epsilon_{t+1}^2 = (y_{t+1} - E[y_{t+1}])^2 = (\\text{values[0]} - 55)^2$  \n",
    "$E[y_{t+3}] = y_{t+2} = $ values[1]  \n",
    "$\\text{VAR}[y_{t+3}] = \\epsilon_{t+2}^2 = (y_{t+2} - E[y_{t+2}])^2 = (\\text{values[1]} - \\text{values[0]})^2$  \n",
    "$E[y_{t+4}] = y_{t+3} = $ values[2]  \n",
    "$\\text{VAR}[y_{t+4}] = \\epsilon_{t+3}^2 = (y_{t+3} - E[y_{t+3}])^2 = (\\text{values[2]} - \\text{values[1]})^2$  \n",
    "\n",
    "<span style='color:red'>**The sims.variance and sims.means values here are not correct beyond the first value. Do not use them.**</span>\n",
    "\n",
    "We can calculate the effective variance of the residuals using the normal_rvs values derived from the seed in the first few cells:  \n",
    "$y_{t+1} = \\sigma_{t+1} z_{t+1} + E[y_{t+1}]$  \n",
    "\n",
    "Therefore\n",
    "\n",
    "$\\sigma_{t+1} = \\frac{(y_{t+1} - E[y_{t+1}])}{z_{t+1}} = \\frac{\\epsilon_{t+1}}{z_{t+1}}$  \n",
    "\n",
    "$\\sigma_{t+1}^2 = (\\frac{\\epsilon_{t+1}}{z_{t+1}})^2$  \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y: [ 0.  1.  1.  2.  3.  5.  8. 13. 21. 34. 55.]\n",
      "-----Means-----\n",
      "means: [55. 55. 55.]\n",
      "y_t+1: 55.0\n",
      "y_t+2: 61.39905867484306\n",
      "y_t+3: 54.74413935810591\n",
      "-----Residuals and variances-----\n",
      "values: [61.399 54.744 59.738]\n",
      "residuals: [ 6.399 -6.655  4.994]\n",
      "eps_t+1: 6.399058674843062\n",
      "eps_t+2: -6.65491931673715\n",
      "eps_t+3: 4.994192159240882\n",
      "vars: [441.    481.948 526.236]\n",
      "residuals**2: [40.948 44.288 24.942]\n",
      "vars based on residuals / normal_rvs: [441.     40.948  44.288]\n",
      "(y_t - E[y_t])^2: 441.0\n",
      "(y_t+1 - E[y_t+1])^2: 40.94795192408424\n",
      "(y_t+2 - E[y_t+2])^2: 44.287951112281256\n"
     ]
    }
   ],
   "source": [
    "np.set_printoptions(linewidth=100, suppress=True, precision=3)\n",
    "test_model = arch.arch_model(y=temp_df['y'], mean='ARX', vol='GARCH', lags=10, p=10, q=0, rescale=False)\n",
    "test_model.distribution = ARCHNormal(seed=seed)\n",
    "test_res = test_model.fix(param_res.params)\n",
    "sims = test_res.forecast(horizon=3, method='simulation', simulations=1, x=None, align='origin')\n",
    "values = sims.simulations.values.squeeze()\n",
    "means = sims.mean.values.squeeze()\n",
    "vars = sims.variance.values.squeeze()\n",
    "residuals = sims._sim.residuals.squeeze()\n",
    "print(f'y: {y}')\n",
    "print('-----Means-----')\n",
    "print(f'means: {means}')\n",
    "print(f'y_t+1: {y[-1]}')\n",
    "print(f'y_t+2: {values[0]}')\n",
    "print(f'y_t+3: {values[1]}')\n",
    "print('-----Residuals and variances-----')\n",
    "print(f'values: {values}')\n",
    "print(f'residuals: {residuals}')\n",
    "print(f'eps_t+1: {values[0] - 55}')\n",
    "print(f'eps_t+2: {values[1] - values[0]}')\n",
    "print(f'eps_t+3: {values[2] - values[1]}')\n",
    "print(f'vars: {vars}')\n",
    "print(f'residuals**2: {residuals**2}')\n",
    "print(f'vars based on residuals / normal_rvs: {(residuals / normal_rvs[:len(residuals)])**2}')\n",
    "print(f'(y_t - E[y_t])^2: {(y[-1] - y[-2])**2}')\n",
    "print(f'(y_t+1 - E[y_t+1])^2: {(values[0] - 55)**2}')\n",
    "print(f'(y_t+2 - E[y_t+2])^2: {(values[1] - values[0])**2}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding exogenuous variable to mean equation\n",
    "\n",
    "lag variable controls number of autoregressive y terms to predict mean while p is the number of residual lag terms to predict residual variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration:      1,   Func. Count:      7,   Neg. LLF: 91.61411266529774\n",
      "Iteration:      2,   Func. Count:     19,   Neg. LLF: 50.75017327614292\n",
      "Iteration:      3,   Func. Count:     28,   Neg. LLF: 36.44096870845732\n",
      "Iteration:      4,   Func. Count:     37,   Neg. LLF: 365.57387348210733\n",
      "Iteration:      5,   Func. Count:     44,   Neg. LLF: 6.151970278719317\n",
      "Iteration:      6,   Func. Count:     51,   Neg. LLF: 6.43027931491611\n",
      "Iteration:      7,   Func. Count:     58,   Neg. LLF: 6.905866135869735\n",
      "Iteration:      8,   Func. Count:     65,   Neg. LLF: 7.741695093707387\n",
      "Iteration:      9,   Func. Count:     72,   Neg. LLF: 1322.8599913432877\n",
      "Iteration:     10,   Func. Count:     83,   Neg. LLF: 2040.7559862722892\n",
      "Iteration:     11,   Func. Count:     93,   Neg. LLF: 240.66942940093227\n",
      "Iteration:     12,   Func. Count:    101,   Neg. LLF: -3.4950800976081555\n",
      "Iteration:     13,   Func. Count:    107,   Neg. LLF: 0.3790684524240998\n",
      "Iteration:     14,   Func. Count:    118,   Neg. LLF: 4298883.246435112\n",
      "Iteration:     15,   Func. Count:    125,   Neg. LLF: 4772291.355366281\n",
      "Iteration:     16,   Func. Count:    132,   Neg. LLF: -2.5045097580834392\n",
      "Iteration:     17,   Func. Count:    139,   Neg. LLF: -4.513087249166446\n",
      "Iteration:     18,   Func. Count:    145,   Neg. LLF: -4.555164294341218\n",
      "Iteration:     19,   Func. Count:    151,   Neg. LLF: 114.09503698933543\n",
      "Iteration:     20,   Func. Count:    158,   Neg. LLF: -4.598099011733307\n",
      "Iteration:     21,   Func. Count:    164,   Neg. LLF: -4.619149385922204\n",
      "Iteration:     22,   Func. Count:    170,   Neg. LLF: -4.626664316484609\n",
      "Iteration:     23,   Func. Count:    176,   Neg. LLF: -4.626687515933295\n",
      "Iteration:     24,   Func. Count:    182,   Neg. LLF: -4.6266882900401525\n",
      "Optimization terminated successfully    (Exit mode 0)\n",
      "            Current function value: -4.6266882900401525\n",
      "            Iterations: 25\n",
      "            Function evaluations: 182\n",
      "            Gradient evaluations: 24\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>AR-X - ARCH Model Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>          <td>y</td>         <th>  R-squared:         </th>  <td>   0.999</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Mean Model:</th>           <td>AR-X</td>        <th>  Adj. R-squared:    </th>  <td>   0.999</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Vol Model:</th>            <td>ARCH</td>        <th>  Log-Likelihood:    </th> <td>   4.62669</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Distribution:</th>        <td>Normal</td>       <th>  AIC:               </th> <td>  0.746623</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>        <td>Maximum Likelihood</td> <th>  BIC:               </th> <td>   2.25955</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th></th>                        <td></td>          <th>  No. Observations:  </th>     <td>10</td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>           <td>Thu, Oct 31 2024</td>  <th>  Df Residuals:      </th>      <td>7</td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>               <td>07:10:31</td>      <th>  Df Model:          </th>      <td>3</td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<caption>Mean Model</caption>\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>      <th>std err</th>      <th>t</th>       <th>P>|t|</th>      <th>95.0% Conf. Int.</th>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Const</th> <td>    0.0685</td>  <td>4.833e-02</td> <td>    1.417</td> <td>    0.156</td>  <td>[-2.624e-02,  0.163]</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>y[1]</th>  <td>    1.6187</td>  <td>5.025e-04</td> <td> 3221.476</td>   <td>0.000</td>      <td>[  1.618,  1.620]</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x</th>     <td>-9.3301e-03</td> <td>6.372e-03</td> <td>   -1.464</td> <td>    0.143</td> <td>[-2.182e-02,3.159e-03]</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<caption>Volatility Model</caption>\n",
       "<tr>\n",
       "      <td></td>        <th>coef</th>     <th>std err</th>      <th>t</th>       <th>P>|t|</th>      <th>95.0% Conf. Int.</th>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>omega</th>    <td>1.4182e-09</td> <td>3.217e-04</td> <td>4.408e-06</td> <td>    1.000</td> <td>[-6.306e-04,6.306e-04]</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>alpha[1]</th> <td>    0.8825</td> <td>    0.465</td> <td>    1.897</td> <td>5.787e-02</td>  <td>[-2.945e-02,  1.794]</td> \n",
       "</tr>\n",
       "</table><br/><br/>Covariance estimator: robust"
      ],
      "text/latex": [
       "\\begin{center}\n",
       "\\begin{tabular}{lclc}\n",
       "\\toprule\n",
       "\\textbf{Dep. Variable:} &         y          & \\textbf{  R-squared:         } &     0.999   \\\\\n",
       "\\textbf{Mean Model:}    &        AR-X        & \\textbf{  Adj. R-squared:    } &     0.999   \\\\\n",
       "\\textbf{Vol Model:}     &        ARCH        & \\textbf{  Log-Likelihood:    } &    4.62669  \\\\\n",
       "\\textbf{Distribution:}  &       Normal       & \\textbf{  AIC:               } &   0.746623  \\\\\n",
       "\\textbf{Method:}        & Maximum Likelihood & \\textbf{  BIC:               } &    2.25955  \\\\\n",
       "\\textbf{}               &                    & \\textbf{  No. Observations:  } &     10      \\\\\n",
       "\\textbf{Date:}          &  Thu, Oct 31 2024  & \\textbf{  Df Residuals:      } &     7       \\\\\n",
       "\\textbf{Time:}          &      07:10:31      & \\textbf{  Df Model:          } &     3       \\\\\n",
       "\\bottomrule\n",
       "\\end{tabular}\n",
       "\\begin{tabular}{lccccc}\n",
       "               & \\textbf{coef} & \\textbf{std err} & \\textbf{t} & \\textbf{P$> |$t$|$} & \\textbf{95.0\\% Conf. Int.}  \\\\\n",
       "\\midrule\n",
       "\\textbf{Const} &       0.0685  &    4.833e-02     &     1.417  &          0.156       &    [-2.624e-02,  0.163]     \\\\\n",
       "\\textbf{y[1]}  &       1.6187  &    5.025e-04     &  3221.476  &        0.000         &     [  1.618,  1.620]       \\\\\n",
       "\\textbf{x}     &  -9.3301e-03  &    6.372e-03     &    -1.464  &          0.143       &   [-2.182e-02,3.159e-03]    \\\\\n",
       "                  & \\textbf{coef} & \\textbf{std err} & \\textbf{t} & \\textbf{P$> |$t$|$} & \\textbf{95.0\\% Conf. Int.}  \\\\\n",
       "\\midrule\n",
       "\\textbf{omega}    &   1.4182e-09  &    3.217e-04     & 4.408e-06  &          1.000       &   [-6.306e-04,6.306e-04]    \\\\\n",
       "\\textbf{alpha[1]} &       0.8825  &        0.465     &     1.897  &      5.787e-02       &    [-2.945e-02,  1.794]     \\\\\n",
       "\\bottomrule\n",
       "\\end{tabular}\n",
       "%\\caption{AR-X - ARCH Model Results}\n",
       "\\end{center}\n",
       "\n",
       "Covariance estimator: robust"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                          AR-X - ARCH Model Results                           \n",
       "==============================================================================\n",
       "Dep. Variable:                      y   R-squared:                       0.999\n",
       "Mean Model:                      AR-X   Adj. R-squared:                  0.999\n",
       "Vol Model:                       ARCH   Log-Likelihood:                4.62669\n",
       "Distribution:                  Normal   AIC:                          0.746623\n",
       "Method:            Maximum Likelihood   BIC:                           2.25955\n",
       "                                        No. Observations:                   10\n",
       "Date:                Thu, Oct 31 2024   Df Residuals:                        7\n",
       "Time:                        07:10:31   Df Model:                            3\n",
       "                                  Mean Model                                  \n",
       "==============================================================================\n",
       "                  coef    std err          t      P>|t|       95.0% Conf. Int.\n",
       "------------------------------------------------------------------------------\n",
       "Const           0.0685  4.833e-02      1.417      0.156   [-2.624e-02,  0.163]\n",
       "y[1]            1.6187  5.025e-04   3221.476      0.000      [  1.618,  1.620]\n",
       "x          -9.3301e-03  6.372e-03     -1.464      0.143 [-2.182e-02,3.159e-03]\n",
       "                               Volatility Model                              \n",
       "=============================================================================\n",
       "                 coef    std err          t      P>|t|       95.0% Conf. Int.\n",
       "-----------------------------------------------------------------------------\n",
       "omega      1.4182e-09  3.217e-04  4.408e-06      1.000 [-6.306e-04,6.306e-04]\n",
       "alpha[1]       0.8825      0.465      1.897  5.787e-02   [-2.945e-02,  1.794]\n",
       "=============================================================================\n",
       "\n",
       "Covariance estimator: robust\n",
       "\"\"\""
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# only one lag on residual with vol='ARCH'\n",
    "test_model = arch.arch_model(y=test_df['y'], x=test_df['x'], mean='ARX', vol='GARCH', lags=1, p=1, q=0, rescale=False)\n",
    "res = test_model.fit()\n",
    "res.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 8\n",
    "start=0 results in first variance of residuals to be nan as mean cannot be computed with missing autoregressive y value when lag>0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [  nan    2.    2.    8.   18.   50.  128.  338.  882. 2312. 6050.]\n"
     ]
    }
   ],
   "source": [
    "y_coeff = 0 # auto-regressive coefficient for mean equation\n",
    "x_coeff = 0 # coefficient for exogenous variable for mean equation\n",
    "res_coeff = 2 # coefficient for residual variance equation\n",
    "test_model.distribution = ARCHNormal(seed=seed)\n",
    "test_res = test_model.fix([0,y_coeff,x_coeff,0,res_coeff])\n",
    "sims = test_res.forecast(horizon=1, start=0, method='simulation', simulations=1, x=x[:,np.newaxis], align='origin')\n",
    "means = sims.mean.values.squeeze()\n",
    "vars = sims.variance.values.squeeze()\n",
    "print(means, vars)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 9\n",
    "start=1 results in first mean and variance of residuals to be computed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y: [ 0.  1.  1.  2.  3.  5.  8. 13. 21. 34. 55.]\n",
      "means: [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "vars: [   2.    2.    8.   18.   50.  128.  338.  882. 2312. 6050.]\n",
      "res_coeff*y[1:]**2: [   2.    2.    8.   18.   50.  128.  338.  882. 2312. 6050.]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_coeff = 0 # auto-regressive coefficient for mean equation\n",
    "x_coeff = 0 # coefficient for exogenous variable for mean equation\n",
    "res_coeff = 2 # coefficient for residual variance equation\n",
    "test_model.distribution = ARCHNormal(seed=seed)\n",
    "test_res = test_model.fix([0,y_coeff,x_coeff,0,res_coeff])\n",
    "sims = test_res.forecast(horizon=1, start=1, method='simulation', simulations=1, x=x[:,np.newaxis], align='origin')\n",
    "means = sims.mean.values.squeeze()\n",
    "vars = sims.variance.values.squeeze()\n",
    "print(f'y: {y}')\n",
    "print(f'means: {means}')\n",
    "print(f'vars: {vars}')\n",
    "print(f'res_coeff*y[1:]**2: {res_coeff*y[1:]**2}')\n",
    "np.allclose(vars, res_coeff*y[1:]**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 10\n",
    "\n",
    "**Testing the alignment of x and y variables when fitting vs forecasting with start not None which produces forecasts 1 time step ahead of each y value**\n",
    "\n",
    "means[i] = $E[\\hat{y}_{t+1}] = \\phi y_{t} + \\beta x_{t}$ so that we get true out of sample forecasted means where $x$ at time $t$ is used  \n",
    "\n",
    "However, for the variance of the out of sample forecasted values e.g. $\\hat{y}_{t+1}$, the in sample residuals are used to compute the variance and the in sample residuals are calculated as $y_t - (y_{t-1} + x_{t})$    \n",
    "\n",
    "So vars[i] = VAR[$\\hat{y}_{t+1}$] $= \\alpha \\epsilon_{t}^2 = \\alpha (y_{t} - E[y_{t}])^2 = \\alpha (y_{t} - (y_{t-1} + x_{t}))^2$\n",
    "\n",
    "In other words, the $x_t$ is used both for the in sample residual at time $t$ and also for the out of sample forecasted mean at time $t+1$ if we do not use the shift in the $x$ variable\n",
    "\n",
    "See [github issue](https://github.com/bashtage/arch/issues/749) and [docs](https://bashtage.github.io/arch/devel//univariate/univariate_forecasting_with_exogenous_variables.html#Conditional-Mean-Alignment-vs.-Forecast-Alignment)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y:                          [0. 7. 2. 5. 9. 3. 8.] with length 7\n",
      "x:                          [1. 2. 4. 2. 1. 9. 6.] with length 7\n",
      "\n",
      "-----Expected means and residuals for fitting-----\n",
      "\n",
      "y[:-1] + x[1:]:             [ 2. 11.  4.  6. 18.  9.] with length 6\n",
      "y[1:] - (y[:-1] + x[1:]):   [  5.  -9.   1.   3. -15.  -1.] with length 6\n",
      "test_res._resid:            [  5.  -9.   1.   3. -15.  -1.] with length 6\n",
      "test_res._resid are in sample residuals calculated as y_t - (y_{t-1} + x_{t})\n",
      "\n",
      "-----Expected means for forecasting-----\n",
      "\n",
      "y + x:                                      [ 1.  9.  6.  7. 10. 12. 14.] with length 7\n",
      "sims.mean:                                  [ 1.  9.  6.  7. 10. 12. 14.]\n",
      "sims.simulations.values - sims.mean:        [    nan   1.524  -9.36    0.75    2.822 -29.266  -1.302]\n",
      "sims._sim.residuals:                        [    nan   1.524  -9.36    0.75    2.822 -29.266  -1.302]\n",
      "sims._sim.residuals are out of sample residuals\n",
      "\n",
      "-----Expected variances-----\n",
      "\n",
      "test_res._resid^2:              [ 25.  81.   1.   9. 225.   1.]\n",
      "sims._sim.residuals/normal_rvs: [ 25.  81.   1.   9. 225.   1.]\n",
      "sims.variance:                  [ nan  25.  81.   1.   9. 225.   1.]\n",
      "The variance of the out of sample forecast are based on the in sample residuals [y_t - (y_{t-1} + x_{t})]^2\n"
     ]
    }
   ],
   "source": [
    "y_coeff = 1.0\n",
    "x_coeff = 1.0\n",
    "res_coeff = 1.0\n",
    "start = 0\n",
    "\n",
    "test_df = pd.DataFrame({'y':np.array([0,7,2,5,9,3,8]).astype(float), 'x':np.array([1,2,4,2,1,9,6]).astype(float)})\n",
    "y = test_df['y'].values\n",
    "x = test_df['x'].values\n",
    "\n",
    "test_model = arch.arch_model(y=test_df['y'], x=test_df['x'], mean='ARX', vol='GARCH', lags=1, p=1, q=0, rescale=False)\n",
    "test_model.distribution = ARCHNormal(seed=seed)\n",
    "test_res = test_model.fix([0., y_coeff, x_coeff, 0., res_coeff])\n",
    "sims = test_res.forecast(horizon=1, start=start, method='simulation', simulations=1, x=x[:,np.newaxis], align='origin')\n",
    "\n",
    "cal_means = y_coeff*y[:-1] + x_coeff*x[1:]\n",
    "eps = y[1:] - cal_means\n",
    "sim_eps = values - means\n",
    "\n",
    "print(f'y:                          {y[start:]} with length {len(y[start:])}')\n",
    "print(f'x:                          {x[start:]} with length {len(x[start:])}')\n",
    "\n",
    "print('\\n-----Expected means and residuals for fitting-----\\n')\n",
    "print(f'y[:-1] + x[1:]:             {cal_means} with length {len(cal_means)}')\n",
    "print(f'y[1:] - (y[:-1] + x[1:]):   {eps} with length {len(eps)}')\n",
    "print(f'test_res._resid:            {test_res._resid} with length {len(test_res._resid)}')\n",
    "print('test_res._resid are in sample residuals calculated as y_t - (y_{t-1} + x_{t})')\n",
    "\n",
    "print('\\n-----Expected means for forecasting-----\\n')\n",
    "print(f'y + x:                                      {y + x} with length {len(y + x)}')\n",
    "print(f'sims.mean:                                  {sims.mean.values.squeeze()}')\n",
    "print(f'sims.simulations.values - sims.mean:        {sims.simulations.values.squeeze() - sims.mean.values.squeeze()}')\n",
    "print(f'sims._sim.residuals:                        {sims._sim.residuals.squeeze()}')\n",
    "print(f'sims._sim.residuals are out of sample residuals')\n",
    "\n",
    "print('\\n-----Expected variances-----\\n')\n",
    "print(f'test_res._resid^2:              {test_res._resid**2}')\n",
    "print(f'sims._sim.residuals/normal_rvs: {(sims._sim.residuals.squeeze()[1:] / normal_rvs[:len(sims._sim.residuals.squeeze())-1])**2}')\n",
    "print(f'sims.variance:                  {sims.variance.values.squeeze()}')\n",
    "print('The variance of the out of sample forecast are based on the in sample residuals [y_t - (y_{t-1} + x_{t})]^2')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 11\n",
    "\n",
    "**Using start=None i.e. forecast starts from the last available observation**\n",
    "\n",
    "<span style='color:red'>**The sims.variance and sims.means values here are not correct beyond the first value. Do not use them.**</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>y</th>\n",
       "      <th>x</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2.0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3.0</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>5.0</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>8.0</td>\n",
       "      <td>6.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>13.0</td>\n",
       "      <td>7.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>21.0</td>\n",
       "      <td>8.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>34.0</td>\n",
       "      <td>9.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>55.0</td>\n",
       "      <td>10.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       y     x\n",
       "0    0.0   0.0\n",
       "1    1.0   1.0\n",
       "2    1.0   2.0\n",
       "3    2.0   3.0\n",
       "4    3.0   4.0\n",
       "5    5.0   5.0\n",
       "6    8.0   6.0\n",
       "7   13.0   7.0\n",
       "8   21.0   8.0\n",
       "9   34.0   9.0\n",
       "10  55.0  10.0"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_df = pd.DataFrame({'y':np.array([0,1,1,2,3,5,8,13,21,34,55]).astype(float), 'x':np.arange(11).astype(float)})\n",
    "y = test_df['y'].values\n",
    "x = test_df.x.values\n",
    "test_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration:      1,   Func. Count:      7,   Neg. LLF: 91.61411266529774\n",
      "Iteration:      2,   Func. Count:     19,   Neg. LLF: 50.75017327614292\n",
      "Iteration:      3,   Func. Count:     28,   Neg. LLF: 36.44096870845732\n",
      "Iteration:      4,   Func. Count:     37,   Neg. LLF: 365.57387348210733\n",
      "Iteration:      5,   Func. Count:     44,   Neg. LLF: 6.151970278719317\n",
      "Iteration:      6,   Func. Count:     51,   Neg. LLF: 6.43027931491611\n",
      "Iteration:      7,   Func. Count:     58,   Neg. LLF: 6.905866135869735\n",
      "Iteration:      8,   Func. Count:     65,   Neg. LLF: 7.741695093707387\n",
      "Iteration:      9,   Func. Count:     72,   Neg. LLF: 1322.8599913432877\n",
      "Iteration:     10,   Func. Count:     83,   Neg. LLF: 2040.7559862722892\n",
      "Iteration:     11,   Func. Count:     93,   Neg. LLF: 240.66942940093227\n",
      "Iteration:     12,   Func. Count:    101,   Neg. LLF: -3.4950800976081555\n",
      "Iteration:     13,   Func. Count:    107,   Neg. LLF: 0.3790684524240998\n",
      "Iteration:     14,   Func. Count:    118,   Neg. LLF: 4298883.246435112\n",
      "Iteration:     15,   Func. Count:    125,   Neg. LLF: 4772291.355366281\n",
      "Iteration:     16,   Func. Count:    132,   Neg. LLF: -2.5045097580834392\n",
      "Iteration:     17,   Func. Count:    139,   Neg. LLF: -4.513087249166446\n",
      "Iteration:     18,   Func. Count:    145,   Neg. LLF: -4.555164294341218\n",
      "Iteration:     19,   Func. Count:    151,   Neg. LLF: 114.09503698933543\n",
      "Iteration:     20,   Func. Count:    158,   Neg. LLF: -4.598099011733307\n",
      "Iteration:     21,   Func. Count:    164,   Neg. LLF: -4.619149385922204\n",
      "Iteration:     22,   Func. Count:    170,   Neg. LLF: -4.626664316484609\n",
      "Iteration:     23,   Func. Count:    176,   Neg. LLF: -4.626687515933295\n",
      "Iteration:     24,   Func. Count:    182,   Neg. LLF: -4.6266882900401525\n",
      "Optimization terminated successfully    (Exit mode 0)\n",
      "            Current function value: -4.6266882900401525\n",
      "            Iterations: 25\n",
      "            Function evaluations: 182\n",
      "            Gradient evaluations: 24\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>AR-X - ARCH Model Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>          <td>y</td>         <th>  R-squared:         </th>  <td>   0.999</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Mean Model:</th>           <td>AR-X</td>        <th>  Adj. R-squared:    </th>  <td>   0.999</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Vol Model:</th>            <td>ARCH</td>        <th>  Log-Likelihood:    </th> <td>   4.62669</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Distribution:</th>        <td>Normal</td>       <th>  AIC:               </th> <td>  0.746623</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>        <td>Maximum Likelihood</td> <th>  BIC:               </th> <td>   2.25955</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th></th>                        <td></td>          <th>  No. Observations:  </th>     <td>10</td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>           <td>Wed, Oct 30 2024</td>  <th>  Df Residuals:      </th>      <td>7</td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>               <td>14:26:31</td>      <th>  Df Model:          </th>      <td>3</td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<caption>Mean Model</caption>\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>      <th>std err</th>      <th>t</th>       <th>P>|t|</th>      <th>95.0% Conf. Int.</th>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Const</th> <td>    0.0685</td>  <td>4.833e-02</td> <td>    1.417</td> <td>    0.156</td>  <td>[-2.624e-02,  0.163]</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>y[1]</th>  <td>    1.6187</td>  <td>5.025e-04</td> <td> 3221.476</td>   <td>0.000</td>      <td>[  1.618,  1.620]</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x</th>     <td>-9.3301e-03</td> <td>6.372e-03</td> <td>   -1.464</td> <td>    0.143</td> <td>[-2.182e-02,3.159e-03]</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<caption>Volatility Model</caption>\n",
       "<tr>\n",
       "      <td></td>        <th>coef</th>     <th>std err</th>      <th>t</th>       <th>P>|t|</th>      <th>95.0% Conf. Int.</th>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>omega</th>    <td>1.4182e-09</td> <td>3.217e-04</td> <td>4.408e-06</td> <td>    1.000</td> <td>[-6.306e-04,6.306e-04]</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>alpha[1]</th> <td>    0.8825</td> <td>    0.465</td> <td>    1.897</td> <td>5.787e-02</td>  <td>[-2.945e-02,  1.794]</td> \n",
       "</tr>\n",
       "</table><br/><br/>Covariance estimator: robust"
      ],
      "text/latex": [
       "\\begin{center}\n",
       "\\begin{tabular}{lclc}\n",
       "\\toprule\n",
       "\\textbf{Dep. Variable:} &         y          & \\textbf{  R-squared:         } &     0.999   \\\\\n",
       "\\textbf{Mean Model:}    &        AR-X        & \\textbf{  Adj. R-squared:    } &     0.999   \\\\\n",
       "\\textbf{Vol Model:}     &        ARCH        & \\textbf{  Log-Likelihood:    } &    4.62669  \\\\\n",
       "\\textbf{Distribution:}  &       Normal       & \\textbf{  AIC:               } &   0.746623  \\\\\n",
       "\\textbf{Method:}        & Maximum Likelihood & \\textbf{  BIC:               } &    2.25955  \\\\\n",
       "\\textbf{}               &                    & \\textbf{  No. Observations:  } &     10      \\\\\n",
       "\\textbf{Date:}          &  Wed, Oct 30 2024  & \\textbf{  Df Residuals:      } &     7       \\\\\n",
       "\\textbf{Time:}          &      14:26:31      & \\textbf{  Df Model:          } &     3       \\\\\n",
       "\\bottomrule\n",
       "\\end{tabular}\n",
       "\\begin{tabular}{lccccc}\n",
       "               & \\textbf{coef} & \\textbf{std err} & \\textbf{t} & \\textbf{P$> |$t$|$} & \\textbf{95.0\\% Conf. Int.}  \\\\\n",
       "\\midrule\n",
       "\\textbf{Const} &       0.0685  &    4.833e-02     &     1.417  &          0.156       &    [-2.624e-02,  0.163]     \\\\\n",
       "\\textbf{y[1]}  &       1.6187  &    5.025e-04     &  3221.476  &        0.000         &     [  1.618,  1.620]       \\\\\n",
       "\\textbf{x}     &  -9.3301e-03  &    6.372e-03     &    -1.464  &          0.143       &   [-2.182e-02,3.159e-03]    \\\\\n",
       "                  & \\textbf{coef} & \\textbf{std err} & \\textbf{t} & \\textbf{P$> |$t$|$} & \\textbf{95.0\\% Conf. Int.}  \\\\\n",
       "\\midrule\n",
       "\\textbf{omega}    &   1.4182e-09  &    3.217e-04     & 4.408e-06  &          1.000       &   [-6.306e-04,6.306e-04]    \\\\\n",
       "\\textbf{alpha[1]} &       0.8825  &        0.465     &     1.897  &      5.787e-02       &    [-2.945e-02,  1.794]     \\\\\n",
       "\\bottomrule\n",
       "\\end{tabular}\n",
       "%\\caption{AR-X - ARCH Model Results}\n",
       "\\end{center}\n",
       "\n",
       "Covariance estimator: robust"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                          AR-X - ARCH Model Results                           \n",
       "==============================================================================\n",
       "Dep. Variable:                      y   R-squared:                       0.999\n",
       "Mean Model:                      AR-X   Adj. R-squared:                  0.999\n",
       "Vol Model:                       ARCH   Log-Likelihood:                4.62669\n",
       "Distribution:                  Normal   AIC:                          0.746623\n",
       "Method:            Maximum Likelihood   BIC:                           2.25955\n",
       "                                        No. Observations:                   10\n",
       "Date:                Wed, Oct 30 2024   Df Residuals:                        7\n",
       "Time:                        14:26:31   Df Model:                            3\n",
       "                                  Mean Model                                  \n",
       "==============================================================================\n",
       "                  coef    std err          t      P>|t|       95.0% Conf. Int.\n",
       "------------------------------------------------------------------------------\n",
       "Const           0.0685  4.833e-02      1.417      0.156   [-2.624e-02,  0.163]\n",
       "y[1]            1.6187  5.025e-04   3221.476      0.000      [  1.618,  1.620]\n",
       "x          -9.3301e-03  6.372e-03     -1.464      0.143 [-2.182e-02,3.159e-03]\n",
       "                               Volatility Model                              \n",
       "=============================================================================\n",
       "                 coef    std err          t      P>|t|       95.0% Conf. Int.\n",
       "-----------------------------------------------------------------------------\n",
       "omega      1.4182e-09  3.217e-04  4.408e-06      1.000 [-6.306e-04,6.306e-04]\n",
       "alpha[1]       0.8825      0.465      1.897  5.787e-02   [-2.945e-02,  1.794]\n",
       "=============================================================================\n",
       "\n",
       "Covariance estimator: robust\n",
       "\"\"\""
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# only one lag on residual with vol='ARCH'\n",
    "test_model = arch.arch_model(y=test_df['y'], x=test_df['x'], mean='ARX', vol='GARCH', lags=1, p=1, q=0, rescale=False)\n",
    "res = test_model.fit()\n",
    "res.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y:                     [ 1.  1.  2.  3.  5.  8. 13. 21. 34. 55.] with length 10\n",
      "x:                     [ 1.  2.  3.  4.  5.  6.  7.  8.  9. 10.] with length 10\n",
      "values:                [31.645 17.504 29.429 12.343  2.887] with length 5\n",
      "-----Means-----\n",
      "means:                 [23.5   21.05  27.315  9.595  4.978] with length 5\n",
      "means from residuals:  [23.5   23.493 26.251 10.229  5.803]\n",
      "E[y_t+1] =             23.5\n",
      "E[y_t+2] =             23.49340160238375\n",
      "E[y_t+3] =             26.251192265196096\n",
      "E[y_t+4] =             10.228843007992765\n",
      "E[y_t+5] =             5.8027966281988785\n",
      "-----Residuals and variances-----\n",
      "residuals:             [ 8.145 -5.989  3.178  2.114 -2.916] with length 5\n",
      "eps_t+1:               8.144672007945832\n",
      "eps_t+2:               -5.989427385063433\n",
      "eps_t+3:               3.1782844281131233\n",
      "eps_t+4:               2.1138124193368313\n",
      "eps_t+5:               -2.9161949242087943\n",
      "vars:                  [714.42   97.466  26.709   7.455   2.905]\n",
      "var from residuals:    [33.168 17.937  5.051  2.234  4.252]\n",
      "residuals/normal_rvs:  [714.42   33.168  17.937   5.051   2.234]\n",
      "(y_t - E[y_t])^2:      714.4199999999998\n",
      "(y_t+1 - E[y_t+1])^2:  33.167841058508195\n",
      "(y_t+2 - E[y_t+2])^2:  17.936620200473897\n",
      "(y_t+3 - E[y_t+3])^2:  5.050745952993181\n",
      "(y_t+4 - E[y_t+4])^2:  2.234101472071314\n",
      "(y_t+5 - E[y_t+5])^2:  4.252096417990568\n"
     ]
    }
   ],
   "source": [
    "y_coeff = 0.3\n",
    "x_coeff = 0.7\n",
    "res_coeff = 0.5\n",
    "\n",
    "x_forecast = [10,20,30,2,3]\n",
    "test_model.distribution = ARCHNormal(seed=seed)\n",
    "test_res = test_model.fix([0., y_coeff, x_coeff, 0., res_coeff])\n",
    "sims = test_res.forecast(horizon=5, method='simulation', simulations=1, x=x_forecast, align='origin')\n",
    "values = sims.simulations.values.squeeze()\n",
    "residuals = sims._sim.residuals.squeeze()\n",
    "means = sims.mean.values.squeeze()\n",
    "vars = sims.variance.values.squeeze()\n",
    "\n",
    "values = sims.simulations.values.squeeze()\n",
    "print(f'y:                     {y[start:]} with length {len(y[start:])}')\n",
    "print(f'x:                     {x[start:]} with length {len(x[start:])}')\n",
    "print(f'values:                {values} with length {len(values)}')\n",
    "print('-----Means-----')\n",
    "print(f'means:                 {means} with length {len(means)}')\n",
    "print(f'means from residuals:  {values - residuals}')\n",
    "print(f'E[y_t+1] =             {y_coeff*y[-1] + x_coeff*x_forecast[0]}')\n",
    "print(f'E[y_t+2] =             {y_coeff*values[0] + x_coeff*x_forecast[1]}')\n",
    "print(f'E[y_t+3] =             {y_coeff*values[1] + x_coeff*x_forecast[2]}')\n",
    "print(f'E[y_t+4] =             {y_coeff*values[2] + x_coeff*x_forecast[3]}')\n",
    "print(f'E[y_t+5] =             {y_coeff*values[3] + x_coeff*x_forecast[4]}')\n",
    "print('-----Residuals and variances-----')\n",
    "print(f'residuals:             {residuals} with length {len(residuals)}')\n",
    "print(f'eps_t+1:               {values[0] - (y_coeff*y[-1] + x_coeff*x_forecast[0])}')\n",
    "print(f'eps_t+2:               {values[1] - (y_coeff*values[0] + x_coeff*x_forecast[1])}')\n",
    "print(f'eps_t+3:               {values[2] - (y_coeff*values[1] + x_coeff*x_forecast[2])}')\n",
    "print(f'eps_t+4:               {values[3] - (y_coeff*values[2] + x_coeff*x_forecast[3])}')\n",
    "print(f'eps_t+5:               {values[4] - (y_coeff*values[3] + x_coeff*x_forecast[4])}')\n",
    "print(f'vars:                  {vars}')\n",
    "print(f'var from residuals:    {res_coeff*residuals**2}')\n",
    "print(f'residuals/normal_rvs:  {(residuals / normal_rvs[:len(residuals)])**2}')\n",
    "print(f'(y_t - E[y_t])^2:      {res_coeff*(y[-1] - (y_coeff*y[-2] + x_coeff*x[-1]))**2}')\n",
    "print(f'(y_t+1 - E[y_t+1])^2:  {res_coeff*(values[0] - (y_coeff*y[-1] + x_coeff*x_forecast[0]))**2}')\n",
    "print(f'(y_t+2 - E[y_t+2])^2:  {res_coeff*(values[1] - (y_coeff*values[0] + x_coeff*x_forecast[1]))**2}')\n",
    "print(f'(y_t+3 - E[y_t+3])^2:  {res_coeff*(values[2] - (y_coeff*values[1] + x_coeff*x_forecast[2]))**2}')\n",
    "print(f'(y_t+4 - E[y_t+4])^2:  {res_coeff*(values[3] - (y_coeff*values[2] + x_coeff*x_forecast[3]))**2}')\n",
    "print(f'(y_t+5 - E[y_t+5])^2:  {res_coeff*(values[4] - (y_coeff*values[3] + x_coeff*x_forecast[4]))**2}')\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "signature",
   "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}