{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c4437c6a-740a-458c-9a03-acf22e7d47c3",
   "metadata": {},
   "source": [
    "1. Raw data manipulation\n",
    "\n",
    "Download the csv file containing the data \"MaxTemperaturesParis\" and place it in a folder. Complete the code to load this csv file and extract the time stamps and measurements as two np.array objects."
   ]
  },
  {
   "cell_type": "raw",
   "id": "c3c88124-bcd2-479f-96ae-db5b5350979d",
   "metadata": {},
   "source": [
    "import csv\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "path_to_data = '???/'\n",
    "file = path_to_data + \"???.csv\"\n",
    "\n",
    "####################################################\n",
    "def readMyFile(filename):\n",
    "    dates = []\n",
    "    scores = ???\n",
    "    with open(filename) as csvDataFile:\n",
    "        csvReader = csv.reader(csvDataFile)\n",
    "        for ???:\n",
    "            dates.append(float(row[0]))\n",
    "            ???\n",
    "    return ???, np.array(scores)\n",
    "####################################################\n",
    "\n",
    "time,meas = readMyFile(???)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d02287dc-a283-420e-92f2-c933ccb80f12",
   "metadata": {},
   "source": [
    "Plot the time trends, superposing the data point in green triangles and the connected data in joined blue lines"
   ]
  },
  {
   "cell_type": "raw",
   "id": "88fc2edc-1c74-4e41-93a7-8a18a0cdbe97",
   "metadata": {},
   "source": [
    "# Fig 2.1\n",
    "plt.plot(time, ???, 'g^', ???, ???, 'b-')\n",
    "plt.ylabel('Temperature (°C)')\n",
    "plt.xlabel('Year')\n",
    "plt.grid(linestyle=':', linewidth=0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c8ce896-d9f7-4789-b3ee-67afefad8dd2",
   "metadata": {},
   "source": [
    "2. Time series decomposition:\n",
    "\n",
    "- choose the length l of the decomposition you want to achieve;\n",
    "- build the associated trajectory (Hankel) matrix Yl;\n",
    "- perform Singular Value Decomposition on this matrix;\n",
    "- plot the singular values and their amplitude;\n",
    "- plot the lead right singular vectors and indicate their respective contributions."
   ]
  },
  {
   "cell_type": "raw",
   "id": "c6ab1202-1075-43db-a307-f78a5541a47c",
   "metadata": {},
   "source": [
    "from scipy.linalg import hankel\n",
    "\n",
    "## Decomposition\n",
    "l = ???\n",
    "N = len(???)\n",
    "Yl = ??? # Build Hankel matrix Y_l\n",
    "U, S, Vh = ??? # Perform SVD on Y_l\n",
    "\n",
    "####################################################\n",
    "def plotSV(S, Vh, mag):\n",
    "    # plot the singular values S and the right singular vectors Vh on the scale +/- mag\n",
    "    l = len(S)-1\n",
    "    plt.suptitle(\"Singular values for SSA\", fontsize=18, y=0.95)\n",
    "    plt.bar(range(1,l+2),???)\n",
    "    plt.ylabel('Amplitude')\n",
    "    plt.xlabel('Orders')\n",
    "    plt.grid(linestyle=':', linewidth=0.5)\n",
    "    plt.show()\n",
    "    plt.figure(figsize=(15, 12))\n",
    "    plt.subplots_adjust(hspace=0.5)\n",
    "    plt.suptitle(\"Lead right singular vectors\", fontsize=18, y=0.95)\n",
    "    for n in range(0,l+1):\n",
    "        # add a new subplot iteratively\n",
    "        ax = plt.subplot(3, 5, n + 1)\n",
    "        ax.plot(???)\n",
    "        ax.set_ylim(-mag,mag)\n",
    "        r_sigma = ??? # Quantifies the contribution of Y_li to Y_l\n",
    "        ax.set_title(\"P.C. %i (%.2f%%)\" % (n+1,r_sigma))\n",
    "        ax.set_xlabel(\"Samples\")\n",
    "####################################################\n",
    "\n",
    "# Fig 2.2 - 2.3\n",
    "plotSV(S, Vh, 0.25)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0abe6fb7-d24a-4838-938f-627c5f847b15",
   "metadata": {},
   "source": [
    "Analyze your results:\n",
    "- how many main components do we have? how do you interpret them?\n",
    "- which amount of the trajectory matrix is contained in these components?\n",
    "- how do you interpret the fact that the 13th singular value is not zero?\n",
    "- how many singular values can be related to the noise spectrum?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cabc5e28-03d0-4c1a-8aca-c48ee95f72a8",
   "metadata": {},
   "source": [
    "3. Reconstruction\n",
    "\n",
    "Fill in the following code to achieve the time series reconstruction of the 1st principal component using the SVD results obtained in the previous section"
   ]
  },
  {
   "cell_type": "raw",
   "id": "fb8ab772-a56e-454b-82fb-27ffeecd2518",
   "metadata": {},
   "source": [
    "# Reconstruction\n",
    "\n",
    "PC = [???] # Selected principal components\n",
    "\n",
    "\n",
    "####################################################\n",
    "def ExtractPC(S,U,Vh,PC):\n",
    "    # Extract the selected principal components PC of the sequence using\n",
    "    #  its singular values decomposition U * S * Vh, in the output matrix Yl_tilde\n",
    "    Yl_tilde = []\n",
    "    iteration = 0\n",
    "    for n in PC:\n",
    "        i = n-1\n",
    "        if iteration == 0:\n",
    "            Yl_tilde = ???\n",
    "            iteration = 1\n",
    "        else:\n",
    "            ???\n",
    "    return Yl_tilde\n",
    "####################################################\n",
    "\n",
    "####################################################\n",
    "# Diagonal averaging\n",
    "def DiagAveraging(Yl_tilde):\n",
    "    # Diagonal averaging method to Hankelize a matrix \"Yl_tilde\" and return the data vector\n",
    "    # N. Golyandina, A. Zhigljavsky, Singular Spectrum Analysis for Time Series, \n",
    "    #  SpringerBriefs in Statistics, Springer Berlin, Heidelberg, 2013.\n",
    "    #  https://link.springer.com/book/10.1007/978-3-642-34913-3\n",
    "    # Code by E. Witrant, May 28th, 2025\n",
    "    L, K = Yl_tilde.shape\n",
    "    L_star = min(L,K)\n",
    "    K_star = max(L,K)\n",
    "    N = L+K-1\n",
    "    if L<K:\n",
    "        Y_star = Yl_tilde\n",
    "    else:\n",
    "        Y_star = np.transpose(Yl_tilde)\n",
    "    y = []\n",
    "    for k in range(1,N+1):\n",
    "        if (k>=1 and k<L_star):\n",
    "            val = 0\n",
    "            for m in range(0,k):\n",
    "                val += Y_star[m][k-1-m]\n",
    "            val = val/k\n",
    "        elif (k>=L_star and k<=K_star):\n",
    "            val = 0\n",
    "            for m in range(0,L_star):\n",
    "                val += Y_star[m][k-1-m]\n",
    "            val = val/L_star\n",
    "        else:\n",
    "            val = 0\n",
    "            for m in range(k-K_star,N-K_star+1):\n",
    "                val += Y_star[m][k-1-m]\n",
    "            val = val/(N-k+1)\n",
    "        y.append(val)\n",
    "    return np.reshape(np.array(y),(len(y),1))\n",
    "####################################################\n",
    "\n",
    "y = DiagAveraging(???)\n",
    "\n",
    "####################################################\n",
    "def plotTrendSeasonality(time, trend, seasonality):\n",
    "    # Plot the trend and seasonality on 2 figures\n",
    "    plt.figure(figsize=(15, 12))\n",
    "    plt.subplots_adjust(hspace=0.5)\n",
    "    plt.suptitle(\"Trend and seasonality\", fontsize=18, y=0.95)\n",
    "    \n",
    "    ax = plt.subplot(2,1,1)\n",
    "    ax.plot(???, ???, 'b-^')\n",
    "    ax.set_title(\"SSA trend time series\")\n",
    "    ax.set_xlabel(\"Year\")\n",
    "    ax.set_ylabel(\"Magnitude\")\n",
    "    ax.grid(linestyle=':', linewidth=0.5)\n",
    "    \n",
    "    ax = plt.subplot(2,1,2)\n",
    "    ax.plot(time, ???, 'b-^')\n",
    "    ax.set_title(\"SSA seasonality time series\")\n",
    "    ax.set_xlabel(\"Year\")\n",
    "    ax.set_ylabel(\"Magnitude\")\n",
    "    ax.grid(linestyle=':', linewidth=0.5)\n",
    "####################################################\n",
    "\n",
    "# Fig 2.4\n",
    "plotTrendSeasonality(time, ???, ???)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39eab40d-a9e2-43da-9a93-f1f4e60ba83b",
   "metadata": {},
   "source": [
    "Group the first three components of the Paris temperature time series and compute \\tilde{y}^{rehan}_t. Compare this result with y_t and plot the residuals."
   ]
  },
  {
   "cell_type": "raw",
   "id": "10cec0e2-4363-43e1-9b17-92dc37e6f90e",
   "metadata": {},
   "source": [
    "PC = ???\n",
    "y_SSA = ???\n",
    "\n",
    "####################################################\n",
    "def plotRawEstimNoise(time,meas,y_SSA):\n",
    "    # Compare the measurement with the estimate y_SSA on a first plot, and show the residuals next.\n",
    "    plt.figure(figsize=(15, 12))\n",
    "    plt.subplots_adjust(hspace=0.5)\n",
    "    plt.suptitle(\"Initial vs. estimated time-series\", fontsize=18, y=0.95)\n",
    "    \n",
    "    ax = plt.subplot(2,1,1)\n",
    "    ax.plot(???, ???, 'k-^', ???, ???, 'r-^')\n",
    "    #    ax.set_ylim(-0.25,0.25)\n",
    "    ax.set_title(\"Raw time series model learning\")\n",
    "    ax.set_xlabel(\"Year\")\n",
    "    ax.set_ylabel(\"Magnitude\")\n",
    "    ax.grid(linestyle=':', linewidth=0.5)\n",
    "    \n",
    "    ax = plt.subplot(2,1,2)\n",
    "    ax.plot(???, ???, 'b-^')\n",
    "    ax.set_title(\"Residuals after model learning\")\n",
    "    ax.set_xlabel(\"Year\")\n",
    "    ax.set_ylabel(\"Magnitude\")\n",
    "    ax.grid(linestyle=':', linewidth=0.5)\n",
    "####################################################\n",
    "\n",
    "# Fig 2.4\n",
    "plotRawEstimNoise(???)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "855838a9-4a38-4db6-bc0a-bf8754fa941f",
   "metadata": {},
   "source": [
    "What can you conclude about this result?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52617b12-0322-4fac-beb3-050e3f18c6cf",
   "metadata": {},
   "source": [
    "Repeat the whole SSA to study the monthly mean carbon dioxide concentration measured at Mauna Loa Observatory, Hawaii, USA, between 1959 and 1991. The data is given in the file \"HawaiiCarbonDioxide\""
   ]
  },
  {
   "cell_type": "raw",
   "id": "b11ca10e-07af-447d-90c1-f8f1d992c24a",
   "metadata": {},
   "source": [
    "file = path_to_data + ???\n",
    "time,meas = readMyFile(file)\n",
    "\n",
    "# Fig 2.6\n",
    "plt.figure(figsize=(15, 10))\n",
    "plt.plot(time, meas, 'k-^')\n",
    "plt.ylabel('CO2 emission (ppm)')\n",
    "plt.xlabel('Year')\n",
    "plt.grid(linestyle=':', linewidth=0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84f2ad69-b5c0-45ca-8b94-987085a20a90",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}