TP3_Time_Series_Kervella.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "TP3_Time_Series_Kervella.ipynb",
      "provenance": [],
      "authorship_tag": "ABX9TyMhM073XRfMfMmwFccU57NT",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/SushiFou/Time-Series-Financial-Data/blob/main/TP3_Time_Series_Kervella.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Eh5vohMvm1Xb"
      },
      "source": [
        "# Time Series for Financial Data - TP n° 3 (Spectral Analysis)\n",
        "---\n",
        "\n",
        "Yann Kervella"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "YUzVZBBmmv-u",
        "outputId": "f5451d2a-34c8-4cf8-b94f-aa9018a71b40"
      },
      "source": [
        "!pip install pyreadr"
      ],
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Collecting pyreadr\n",
            "\u001b[?25l  Downloading https://files.pythonhosted.org/packages/11/ae/74e99f7fe3652f5976acc35543fdf17abe9b887478218c9e277b509c0a45/pyreadr-0.4.0-cp37-cp37m-manylinux2014_x86_64.whl (410kB)\n",
            "\r\u001b[K     |▉                               | 10kB 12.2MB/s eta 0:00:01\r\u001b[K     |█▋                              | 20kB 16.8MB/s eta 0:00:01\r\u001b[K     |██▍                             | 30kB 10.2MB/s eta 0:00:01\r\u001b[K     |███▏                            | 40kB 9.7MB/s eta 0:00:01\r\u001b[K     |████                            | 51kB 4.5MB/s eta 0:00:01\r\u001b[K     |████▉                           | 61kB 4.9MB/s eta 0:00:01\r\u001b[K     |█████▋                          | 71kB 5.1MB/s eta 0:00:01\r\u001b[K     |██████▍                         | 81kB 5.5MB/s eta 0:00:01\r\u001b[K     |███████▏                        | 92kB 5.4MB/s eta 0:00:01\r\u001b[K     |████████                        | 102kB 4.2MB/s eta 0:00:01\r\u001b[K     |████████▉                       | 112kB 4.2MB/s eta 0:00:01\r\u001b[K     |█████████▋                      | 122kB 4.2MB/s eta 0:00:01\r\u001b[K     |██████████▍                     | 133kB 4.2MB/s eta 0:00:01\r\u001b[K     |███████████▏                    | 143kB 4.2MB/s eta 0:00:01\r\u001b[K     |████████████                    | 153kB 4.2MB/s eta 0:00:01\r\u001b[K     |████████████▊                   | 163kB 4.2MB/s eta 0:00:01\r\u001b[K     |█████████████▋                  | 174kB 4.2MB/s eta 0:00:01\r\u001b[K     |██████████████▍                 | 184kB 4.2MB/s eta 0:00:01\r\u001b[K     |███████████████▏                | 194kB 4.2MB/s eta 0:00:01\r\u001b[K     |████████████████                | 204kB 4.2MB/s eta 0:00:01\r\u001b[K     |████████████████▊               | 215kB 4.2MB/s eta 0:00:01\r\u001b[K     |█████████████████▋              | 225kB 4.2MB/s eta 0:00:01\r\u001b[K     |██████████████████▍             | 235kB 4.2MB/s eta 0:00:01\r\u001b[K     |███████████████████▏            | 245kB 4.2MB/s eta 0:00:01\r\u001b[K     |████████████████████            | 256kB 4.2MB/s eta 0:00:01\r\u001b[K     |████████████████████▊           | 266kB 4.2MB/s eta 0:00:01\r\u001b[K     |█████████████████████▌          | 276kB 4.2MB/s eta 0:00:01\r\u001b[K     |██████████████████████▍         | 286kB 4.2MB/s eta 0:00:01\r\u001b[K     |███████████████████████▏        | 296kB 4.2MB/s eta 0:00:01\r\u001b[K     |████████████████████████        | 307kB 4.2MB/s eta 0:00:01\r\u001b[K     |████████████████████████▊       | 317kB 4.2MB/s eta 0:00:01\r\u001b[K     |█████████████████████████▌      | 327kB 4.2MB/s eta 0:00:01\r\u001b[K     |██████████████████████████▍     | 337kB 4.2MB/s eta 0:00:01\r\u001b[K     |███████████████████████████▏    | 348kB 4.2MB/s eta 0:00:01\r\u001b[K     |████████████████████████████    | 358kB 4.2MB/s eta 0:00:01\r\u001b[K     |████████████████████████████▊   | 368kB 4.2MB/s eta 0:00:01\r\u001b[K     |█████████████████████████████▌  | 378kB 4.2MB/s eta 0:00:01\r\u001b[K     |██████████████████████████████▍ | 389kB 4.2MB/s eta 0:00:01\r\u001b[K     |███████████████████████████████▏| 399kB 4.2MB/s eta 0:00:01\r\u001b[K     |████████████████████████████████| 409kB 4.2MB/s eta 0:00:01\r\u001b[K     |████████████████████████████████| 419kB 4.2MB/s \n",
            "\u001b[?25hRequirement already satisfied: pandas>0.24.0 in /usr/local/lib/python3.7/dist-packages (from pyreadr) (1.1.5)\n",
            "Requirement already satisfied: pytz>=2017.2 in /usr/local/lib/python3.7/dist-packages (from pandas>0.24.0->pyreadr) (2018.9)\n",
            "Requirement already satisfied: python-dateutil>=2.7.3 in /usr/local/lib/python3.7/dist-packages (from pandas>0.24.0->pyreadr) (2.8.1)\n",
            "Requirement already satisfied: numpy>=1.15.4 in /usr/local/lib/python3.7/dist-packages (from pandas>0.24.0->pyreadr) (1.19.5)\n",
            "Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.7/dist-packages (from python-dateutil>=2.7.3->pandas>0.24.0->pyreadr) (1.15.0)\n",
            "Installing collected packages: pyreadr\n",
            "Successfully installed pyreadr-0.4.0\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FUICDbDNmzPr"
      },
      "source": [
        "# Spectral analysis on a bivariate time series #\n",
        "\n",
        "\n",
        "## Raw analysis ##\n",
        "\n",
        "We consider two time series derived fom the *realized volatility* of two indices: FTSE and SP500, denoted in the following by $Y_t(1)$ and $Y_t(2)$. \n",
        "Details about the data can be found here:\n",
        "\n",
        " http://realized.oxford-man.ox.ac.uk\n",
        "\n",
        "We first analyze the joint second-order statistics of the two time series and\n",
        "then propose two bivariate dynamic linear models for them, and compare their\n",
        "predictive performence.\n",
        "\n",
        "\n",
        "**1) Load the data in the file** https://perso.telecom-paristech.fr/roueff/edu/tsfd/data/realizedvolatility.Rdata\n",
        "\n",
        "```{r, eval=TRUE}\n",
        "load(url('https://m2:map658@perso.telecom-paristech.fr/roueff/edu/tsfd/data/realizedvolatility.Rdata'))\n",
        "```\n",
        "\n",
        "Apply the following command to have a quick look on the two time series.\n",
        "\n",
        "```{r, eval=FALSE}\n",
        "op <- par(mfrow=c(2,1))\n",
        "plot(as.POSIXct(volsf$Date),volsf$FTSE,type='l',xlab='Date',ylab='FTSE RV')\n",
        "plot(as.POSIXct(volsf$Date),volsf$SPX,type='l',xlab='Date',ylab='SP500 RV')\n",
        "par(op)\n",
        "```\n",
        "From now on, we will actually work with the log volatility:\n",
        "\n",
        "```{r, eval=FALSE}\n",
        "vollog <- data.frame(FTSE=log(volsf$FTSE),SPX=log(volsf$SPX))\n",
        "vollog <- cbind(volsf[,'Date',drop=FALSE],vollog)\n",
        "\n",
        "op <- par(mfrow=c(2,1))\n",
        "plot(as.POSIXct(vollog$Date),vollog$FTSE,type='l',xlab='Date',ylab='Log FTSE RV')\n",
        "plot(as.POSIXct(vollog$Date),vollog$SPX,type='l',xlab='Date',ylab='Log SP500 RV')\n",
        "par(op)\n",
        "```\n",
        "\n",
        "**Have a look at their autocorrelation functions as well as their cross\n",
        "  autocorrelation function (use *acf()* and *ccf()*) and comment.**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "881rq6pknMaR"
      },
      "source": [
        "**2) Draw the periodograms $I_N^{Y(1)}(\\lambda)$ and $I_N^{Y(2)}(\\lambda)$ from the $N$\n",
        "  obserations of $Y(1)$ and $Y(2)$, for $\\lambda=2\\pi k/N$ with $k=1,\\dots,[N/2]$,\n",
        "  computed using *fft()*. Display the periodograms with a log scale in the $y$ axis\n",
        "  by adding the argument** *log='y'* **in the** *plot()* **command**."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oaEtINhcoCsC"
      },
      "source": [
        "## Smoothed periodogram ##\n",
        "\t\n",
        "Periodograms are not consistent estimators of the spectral densities. To obtain\n",
        "consistent estimators, one needs to smooth them over frequencies.\n",
        "  \n",
        "**3) We must first choose the shape of the smoothing kernel. Here are some examples of kernels, which share the same discrete support $\\{-150,\\dots,0,\\dots,150\\}$ :**\n",
        "\n",
        "```{r, eval=FALSE}\n",
        "op <- par(mfrow=c(3,1))\n",
        "plot(kernel(\"daniell\", c(150)))\n",
        "plot(kernel(\"daniell\", c(100,50)))\n",
        "plot(kernel(\"daniell\", c(70,50,30)))\n",
        "par(op)\n",
        "```\n",
        "\n",
        "One can use such kernels to smooth the raw periodogram previously obtained,\n",
        "that is, an estimate of the spectral density is obtained by averaging the\n",
        "periodogram around each frequency. This can be done as follows :\n",
        "\n",
        "\n",
        "```{r, eval=FALSE}\n",
        "sp <- spectrum(cbind(vollog$FTSE,vollog$SPX),kernel=kernel(\"daniell\", c(70,50,30)))\n",
        "```\n",
        "\n",
        "**Discuss the results obtained with kernels of various support lengths (or *bandwidth*).** "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MLFJ02pFn_Xc"
      },
      "source": [
        "**4)Compare with the *raw* periodograms obtained previously.**\n",
        "[Warning: the spectra computed by the default method in *spectrum()* are not normalized using the standard time series convension but the signal processing one. In the time series literature, one usually defines the spectral density $f$ so that \n",
        "$$ \n",
        "\\gamma(h)=\\int_{-\\pi}^{\\pi}\\mathrm{e}^{\\mathrm{i}\\lambda h}\\;f(\\lambda)\\;\\mathrm{d}\\lambda \n",
        "$$ \n",
        "In the signal processing literature, one usually defines the \\emph{power} spectral density $p$ so that \n",
        "$$ \n",
        "\\gamma(h)=\\int_{-1/2}^{1/2}\\mathrm{e}^{2\\mathrm{i}\\pi\\omega h}\\;p(\\omega)\\;\\mathrm{d}\\omega \n",
        "$$ \n",
        "In other words, we have $p(\\omega)=2\\pi\\;f(2\\pi\\omega)$.  ]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "haubgcDIn7rx"
      },
      "source": [
        "**5) The default plot in the method *spectrum()* displays the estimated spectral densities of each time series in *volsf[,-1]*.\n",
        "However, the object returned from *spectrum()* also contains the estimated coherency:**\n",
        "\n",
        "```{r, eval=FALSE}\n",
        "plot(sp,plot.type = \"coherency\")\n",
        "```\n",
        "\n",
        "**Are the two time series more coherent at low frequencies, high frequencies or mid-frequencies ? How can this be interpreted ?** "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BGme9GPNn4i-"
      },
      "source": [
        "# A dynamic linear model #\n",
        "\n",
        "Let us set $\\mathbf{Y}_t=\\begin{bmatrix}Y_t(1)&Y_t(2)\\end{bmatrix}^T$.  Let\n",
        " $X_t$ be a Gaussian AR(1) process with AR coefficient $\\phi$\n",
        "and innovation variance $\\sigma^2$.  We propose here to model the joint dynamics of the two time\n",
        "series $Y_t(1)$ and $Y_t(2)$ by the following equations:\n",
        "\n",
        "\\begin{align*}\n",
        "\\mathbf{Y}_t&= X_t \\mathbf{a} +\\boldsymbol{\\epsilon}_t\\;,\n",
        "\\end{align*}\n",
        "\n",
        "where $\\mathbf{a}$ is a column vector of the form $\\begin{bmatrix}1&\n",
        "a\\end{bmatrix}^T$ and $(\\boldsymbol{\\epsilon}_t)$ is a Gaussian bivariate white\n",
        "noise with covariance matrix Cov$(\\boldsymbol{\\epsilon}_0)=R$ which is independent of $(X_t)$. \n",
        "\n",
        "**5) Recall the form of the spectral density $f^X$ of $(X_t)$. Express the\n",
        "  spectral densities $f^{Y(1)}$ and $f^{Y(2)}$ of $Y_t(1)$ and $Y_t(2)$ and the\n",
        "  coherency $C^Y$ between them using $f^X$, $a$ and $R$. Code an R function\n",
        "  returning $f^X(\\lambda)$, $f^{Y(1)}(\\lambda)$, $f^{Y(2)}(\\lambda)$, and\n",
        "  $C^Y(\\lambda)$, with inputs : a list of frequencies $\\lambda$'s and the\n",
        "  parameters $\\phi$, $\\sigma$, $a$ and $R$. Plot the graph of $C$ for\n",
        "  $\\phi=.5$, $\\sigma=1$, $a=.5$ and $R$ equal to the identity matrix.**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "inuqOvZOn0d9"
      },
      "source": [
        "**6) How does the previous graph of $C$ evolve as $\\phi$ increases towards 1 ?\n",
        "  What is the behavior of the coherency when $f^X(\\lambda)\\gg R$ ? How does\n",
        "  this make a plausible model for the data at hand ?**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NmfiD5T0nyWc"
      },
      "source": [
        "**7) Show that this model is a *dynamic linear model* (DLM) and provide the\n",
        "  equations defining the dynamics of the model.**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WjWPFUi-nwSj"
      },
      "source": [
        "**8) Code an R function returning a Gaussian sample of this DLM with inputs :\n",
        "  the sample length, and the parameters $\\phi$, $\\sigma$, $a$ and $cR$, where\n",
        "  $cR$ is a matrix such that $cR^T\\,cR=R$. The initial value of the state\n",
        "  variable will be drawn according to the stationary distribution. Based on a\n",
        "  simulated $2^{10}$ bivariate sample, with the parameters $\\phi=.5$, $\\sigma=1$,\n",
        "  $a=.5$ and $R$ equal to the identity matrix, perform a spectral analysis\n",
        "  similar to that performed on the log volatility data.**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FpUPrVeOniVb"
      },
      "source": [
        "# Fitting the model #\n",
        "\n",
        "The previous analysis supports the idea of using the proposed DLM for the *centered* log volatility data.\n",
        "We will relay on the following R package to perform computations related to the Kalman algorithm. \n",
        "\n",
        "```{r, eval=FALSE} \n",
        "require(astsa) \n",
        "``` \n",
        "\n",
        "We will mainly use the command\n",
        "*Kfilter0()* in this package.  The inputs of *Kfilter0()* include the DLM\n",
        "parameters $A$, $\\mu_0$, $\\Sigma_0$, $\\Phi$, $cQ$ and $cR$.  The notational\n",
        "convention for these parameters is given by writing the state equation and the\n",
        "observation equation as \n",
        "\n",
        "\\begin{align*} \n",
        "X_t &= \\Phi\\, X_{t-1} + W_t\\;,\\\\ \n",
        "Y_t &= A\\, X_t + V_t\\;, \n",
        "\\end{align*} \n",
        "where $X_0 \\sim \\mathcal{N}(\\mu_0, \\Sigma_0)$, $W_t$ are\n",
        "iid $\\mathcal{N}(0,cQ^T cQ)$, and $V_t$ are iid $\\mathcal{N}(0,cR^T cR)$.\n",
        "\n",
        "The other inputs are the data *y* and the number of observations *num*\n",
        "(corresponding to the rows of *y*).\n",
        "\n",
        "**9) Code an R function which returns $A$, $\\mu_0$, $\\Sigma_0$, $\\Phi$ and $cQ$ of\n",
        "the DLM of Question 7) from the inputs $\\phi$, $\\sigma$ and $a$. (Note that\n",
        "$cR$ is the same.) Using *Kfilter0()*, code a funtion *llike()*, whose input is\n",
        "a vector *para* that contains the parameters $\\phi$, $\\sigma$, $a$ and $cR$,\n",
        "and which returns the corresponding negated log likelihood. The observations\n",
        "used to compute this likelihood needs to be called through an external\n",
        "variable, say *yest*.**\n",
        "\n",
        "The correspondance between *para* and the parameters $\\phi$, $\\sigma$, $a$ and\n",
        "$cR$ can be done as follows\n",
        "\n",
        "```{r, eval=FALSE} \n",
        "a <- para[1] \n",
        "phi <- para[3] \n",
        "sigma <- para[4] \n",
        "cR <- matrix(para[4:7],nrow=2,ncol=2) \n",
        "```\n",
        "With the function *llike()* at hand, one can use *optime()* to numerically minimize it:\n",
        "\n",
        "```{r, eval=FALSE} \n",
        "est <- optim(init.par, llike, NULL, method=\"BFGS\",\n",
        "   hessian=TRUE, control=list(trace=1,REPORT=1))\n",
        "```\n",
        "\n",
        "Here *optim()* uses the iterative \"BFGS\" algorithm and starts from the initial parameter *init.par*. \n",
        "The output *est* contains the resulting locally minimizing parameter *est$par*."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QEjY2upGnczT"
      },
      "source": [
        "**10) Test the previous method on the simulated data of question 8). Use an increasing number of observations, from a sample size $2^6$ to $2^{10}$. Plot the evolution of the estimated $a$ as a function of the sample size, and do the same for $phi$.** [Comment: iterative algorithms can be sensitive to the initial parameter. When the true parameters of the data is known, it is possible to set *init.par* to the true parameters. However it is interesting to compare the output when starting from different initial parameters.]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eN5_YmgXnZpy"
      },
      "source": [
        "**11) Use the $2^{10}$ first samples of *vollog* to estimate the parameters of the DLM that best fits these data. Before fitting the model the data must be centered:**\n",
        "\n",
        "```{r, eval=FALSE} \n",
        "yest <- vollog[1:2^10,-1]-colMeans(vollog[1:2^10,-1])\n",
        "```\n",
        "\n",
        "**Compare the spectral densities and the coherency of the fitted model with the spectra obtained in Questions 3) and 5).**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "b_IU-0fGnSHW"
      },
      "source": [
        "**12) Use the fitted parameters to compute one step ahead predictors of the log volatilities of FTSE and SPX given their past.**\n",
        "\n",
        "[Comment: the output of *Kfilter0()* contains one step ahead predictors of the state variable, called *xp*. It is then easy to deduce predictors of the observed variable using the matrix $A$.]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BtUsJIMBnOp1"
      },
      "source": [
        "## Tiebraker question ##\n",
        "\n",
        "**13) We derived predictors of the *log* volatility in Question 12). How should\n",
        "we deduce predictors of the volatility ?  Compare the prediction performance of\n",
        "these predictors with predictors obtained *individually* from each\n",
        "time series FTSE and SPX.**"
      ]
    }
  ]
}