introducing-ssa-for-time-series-decomposition.ipynb

{"cells":[{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5"},"cell_type":"markdown","source":"# Decomposing Time Series Data With Singular-Spectrum Analysis\n\nThis tutorial will introduce the technique of **singular-spectrum analysis (SSA)**, a method of time series decomposition that has received little attention on Kaggle. This is surprising, as there is ongoing demand for techniques that can be used to explore, analyse and process time series data, the recent [Recruit Restaurant Visitor Forecasting](https://www.kaggle.com/c/recruit-restaurant-visitor-forecasting) competition being one example.\n\nPut simply, SSA decomposes a time series into a set of summable components that are grouped together and interpreted as *trend*, *periodicity* and *noise*. SSA emphasises **separability** of the underlying components, and can readily separate periodicities that occur on different time scales, even in very noisy time series data. The original time series is recovered by summing together all of its components.\n\nThus, SSA can be used to analyse and reconstruct a time series with or without different components as desired. For example, you could apply SSA to:\n* construct a smoothed version of a time series using a small subset of its components; \n* investigate a time series' periodic components to understand the underlying processes that generated the time series;\n* reconstruct the original time series without its periodic components;\n* remove all trend and periodic components from the series, leaving just the 'noise', which may be meaningful in and of itself...\n\nand so on and so forth. Unlike the commonly used autoregressive integrated moving average (ARIMA) method, SSA makes no assumptions about the nature of the time series, and has just a single adjustable (and easily interpretable) parameter. Further, the base algorithm requires little more than a few lines of `numpy` linear algebra to implement, so you can easily roll your own bespoke SSA implementation if you desire.\n\n## The Purpose of This Kernel\nI intend this kernel to not only introduce and demonstrate the SSA method, but also present its theory and Python implementation. I will present theory and its corresponding code in tandem, bridging the gap between mathematical underpinning and practical algorithmic implementation. While R and Python packages for SSA are available (e.g. [Rssa](https://cran.r-project.org/package=Rssa) and [pySSA](https://github.com/aj-cloete/pySSA)), I feel there is more to be gained when the inner workings of a given statistical/analytical technique are understood, instead of relying on a few subroutines from a package.\n\nSSA will be applied to mock time series data to demonstrate its key abilities, illustrate where it excels and falls down, and document potential pitfalls along the way. SSA will then be applied to a real-world time series taken from the dataset, [MotionSense Dataset : Smartphone Sensor Data](https://www.kaggle.com/malekzadeh/motionsense-dataset). This will not be an extensive analysis of the MotionSense dataset, but will instead highlight the ability of SSA to extract underlying periodicities in a noisy, complex time series.\n\nOverall, I hope that these elementary studies of the SSA technique will provide a solid foundation for others to successfully apply SSA to real-world time series data.\n\n## Outline\nIn [Section 1](#Section1), I will define a simple, illustrative toy time series that consists of trend, periodic and noise elements. This time series will be used to demonstrate the workings of the SSA method over the next several sections.\n\n[Section 2](#Section2) introduces the SSA method and translates its formal steps into a practical Python implementation. As we go, I will apply SSA to our toy time series, and intermediate results will be plotted and discussed. At the heart of this section is the acclaimed singular value decomposition. \n\n[Section 3](#Section3) focuses on the notion of separability in SSA, and how it is quantified. Once again, the formal mathematics will be translated into Python, and applied to the toy time series.\n\n[Section 4](#Section4) will bring together the code developed in Sections 2 and 3, creating a simple Python SSA class that will be used in the rest of the kernel (and future kernels).\n\n[Section 5](#Section5) investigates the effect of the window length parameter on the decomposition of the toy time series. \n\n[Section 6](#Section6) applies SSA to real-world time series data from the [MotionSense](https://www.kaggle.com/malekzadeh/motionsense-dataset) dataset. I also propose a hypothetical use for SSA with this dataset. \n\n[Section 7](#Section7) contains a few concluding words.\n\n## Contents\n[1. A Toy Time Series](#Section1) <br>\n[2. Introducing the SSA Method](#Section2) <br>\n&nbsp; &nbsp; &nbsp;[2.1 From Time Series to Trajectory Matrix](#Section2.1) <br>\n&nbsp; &nbsp; &nbsp;[2.2 Decomposing the Trajectory Matrix](#Section2.2) <br>\n&nbsp; &nbsp; &nbsp; &nbsp;&nbsp;&nbsp; [2.2.1 The $\\mathbf{U}$ Matrix](#Section2.2.1) <br>\n&nbsp; &nbsp; &nbsp; &nbsp;&nbsp;&nbsp; [2.2.2 The $\\mathbf{V}$ Matrix](#Section2.2.2) <br>\n&nbsp; &nbsp; &nbsp; &nbsp;&nbsp;&nbsp; [2.2.3 The $\\mathbf{\\Sigma}$ Matrix](#Section2.2.3) <br>\n&nbsp; &nbsp; &nbsp; &nbsp;&nbsp;&nbsp; [2.2.4 The Rank of the Trajectory Matrix](#Section2.2.4) <br>\n&nbsp; &nbsp; &nbsp; &nbsp;&nbsp;&nbsp; [2.2.5 Putting it all Together](#Section2.2.5) <br>\n&nbsp; &nbsp; &nbsp;[2.3 Reconstructing the Time Series](#Section2.3) <br>\n[3. Time Series Component Separation and Grouping](#Section3) <br>\n[4. A Python Class for SSA](#Section4) <br>\n[5. The Window Length](#Section5) <br>\n&nbsp; &nbsp; &nbsp; [5.1 $L=2$](#Section5.1) <br>\n&nbsp; &nbsp; &nbsp; [5.2 $L=5$](#Section5.2) <br>\n&nbsp; &nbsp; &nbsp; [5.3 $L=20$](#Section5.3) <br>\n&nbsp; &nbsp; &nbsp; [5.4 $L=40$](#Section5.4) <br>\n&nbsp; &nbsp; &nbsp; [5.5 $L=60$](#Section5.4) <br>\n[6. MotionSense: Applying SSA to Accelerometer Data](#Section6) <br>\n&nbsp; &nbsp; &nbsp;[6.1 Loading the Data](#Section6.1) <br>\n&nbsp; &nbsp; &nbsp;[6.2 Plotting the Data](#Section6.2) <br>\n&nbsp; &nbsp; &nbsp;[6.3 Decomposing the Time Series with SSA](#Section6.3) <br>\n&nbsp; &nbsp; &nbsp;[6.4 Using SSA to Extract an Individual's Walking 'Signature'](#Section6.4) <br>\n&nbsp; &nbsp; &nbsp;&nbsp; &nbsp; &nbsp;[6.4.1 Data Extraction and Processing](#Section6.4.1) <br>\n&nbsp; &nbsp; &nbsp;&nbsp; &nbsp; &nbsp;[6.4.2 Comparison of Walk 'Signatures'](#Section6.4.2) <br>\n[7. Some Final Words](#Section7) <br>\n## Python Setup\nThis kernel requires only the standard `numpy`, `matplotlib` and `pandas` modules:"},{"metadata":{"collapsed":true,"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"code","source":"# Load the usual suspects:\nimport numpy as np\nfrom numpy import pi\nimport matplotlib.pyplot as plt\nimport pandas as pd\n\n# Fiddle with figure settings here:\nplt.rcParams['figure.figsize'] = (10,8)\nplt.rcParams['font.size'] = 14\nplt.rcParams['image.cmap'] = 'plasma'\nplt.rcParams['axes.linewidth'] = 2\n# Set the default colour cycle (in case someone changes it...)\nfrom cycler import cycler\ncols = plt.get_cmap('tab10').colors\nplt.rcParams['axes.prop_cycle'] = cycler(color=cols)\n\n# A simple little 2D matrix plotter, excluding x and y labels.\ndef plot_2d(m, title=\"\"):\n    plt.imshow(m)\n    plt.xticks([])\n    plt.yticks([])\n    plt.title(title)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"18a1b79d-68a4-42cc-8b99-002c8a0b0717","_uuid":"efef64e6a277270f5fdf52579e389fe3fa2a5d58"},"cell_type":"markdown","source":"# 1. A Toy Time Series <a name=\"Section1\"></a>\nTo demonstrate the workings and concepts of SSA, we'll first define an arbitrary, toy time series, $F = \\{f_0, f_1, \\ldots, f_{N-1}\\}$, containing trend, periodic and noise components:\n\n$$f_t = 0.001 \\times (t - 100)^2 + 2\\sin(\\frac{2\\pi t}{p_1}) + 0.75\\sin(\\frac{2\\pi t}{p_2}) + \\text{Rand}\\{-1,1\\}$$\n\nwhere $t = \\{0, 1,\\ldots, N-1\\}$ is a discrete time moment, $p_1$ and $p_2$ are set time periods and $\\text{Rand}\\{-1,1\\}$ is a random number uniformly distributed between –1 and 1.\n\nThe $0.001 \\times (t - 100)^2$ term defines the (parabolic) trend of the series,  $2\\sin(\\frac{2\\pi t}{p_1})$ and $0.75\\sin(\\frac{2\\pi t}{p_2})$ are two periodic components with differing periodicities and amplitudes, while $\\text{Rand}\\{-1,1\\}$ introduces noise.\nLet's take a look:"},{"metadata":{"collapsed":true,"_cell_guid":"bbca025b-7baa-4959-b0c5-6c16b545b00e","_uuid":"e1b2b3e64c57d6cb9fb4bbc97abfeee20155b524","trusted":false},"cell_type":"code","source":"N = 200 # The number of time 'moments' in our toy series\nt = np.arange(0,N)\ntrend = 0.001 * (t - 100)**2\np1, p2 = 20, 30\nperiodic1 = 2 * np.sin(2*pi*t/p1)\nperiodic2 = 0.75 * np.sin(2*pi*t/p2)\n\nnp.random.seed(123) # So we generate the same noisy time series every time.\nnoise = 2 * (np.random.rand(N) - 0.5)\nF = trend + periodic1 + periodic2 + noise\n\n# Plot everything\nplt.plot(t, F, lw=2.5)\nplt.plot(t, trend, alpha=0.75)\nplt.plot(t, periodic1, alpha=0.75)\nplt.plot(t, periodic2, alpha=0.75)\nplt.plot(t, noise, alpha=0.5)\nplt.legend([\"Toy Series ($F$)\", \"Trend\", \"Periodic #1\", \"Periodic #2\", \"Noise\"])\nplt.xlabel(\"$t$\")\nplt.ylabel(\"$F(t)$\")\nplt.title(\"The Toy Time Series and its Components\");","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"8acc571a-3891-4b14-9c8b-3f52483e7c9d","_uuid":"c2a3664ae4c118efd7219366ca1426b86ec14d70"},"cell_type":"markdown","source":"Based on visual inspection, the presence of trend (orange), periodic (green) and noise (purple) components in the time series is clear. However, the second periodic component (red) in the toy series is not apparent. Could we use SSA to recover the trend, periodic and noise components from the toy time series?\n\n# 2. Introducing the SSA Method<a name=\"Section2\"></a>\nThe notation in this section will mostly follow that used in the book, [*Analysis of Time Series Structure: SSA and Related Techniques*](https://www.crcpress.com/Analysis-of-Time-Series-Structure-SSA-and-Related-Techniques/Golyandina-Nekrutkin-Zhigljavsky/p/book/9781584881940) by N. Golyandina, V. Nekrutkin and A. Zhigljavsky—an excellent reference containing more proofs and details about SSA than you could possibly ever need to know. Here, we'll go a step further and develop a Python implementation for SSA. For consistency between the mathematical notation and the code, all matrix and summation indexing with begin at 0 instead of 1.\n\n## 2.1 From Time Series to Trajectory Matrix<a name=\"Section2.1\"></a> \nThe first step of SSA is to map the time series $F$ to a sequence of multi-dimensional lagged vectors. Let an integer $L$ be the **window length**, $2 \\le L \\le N/2$. We form a 'window', given by the subseries $\\{f_i, \\ f_{i+1}, \\ldots , \\ f_{i+L-1}\\}$, for  $i=0,\\ldots,N-L$. We slide this window along the time series, forming a column vector, $X_i$, for each subseries. That is, we have\n\n\\begin{align*}\nX_0 & = (f_0, \\ f_1, \\ f_2,  \\ldots, \\ f_{L-1} )^{\\text{T}} \\\\\nX_1 & = (f_1, \\ f_2, \\ f_3,  \\ldots, \\ f_L )^{\\text{T}} \\\\\nX_2 & = (f_2, \\ f_3, \\ f_4,  \\ldots, \\ f_{L+1} )^{\\text{T}} \\\\\nX_3 & = (f_3, \\ f_4, \\ f_5,  \\ldots, \\ f_{L+2} )^{\\text{T}} \\\\\n& \\quad \\quad \\quad  \\vdots \\\\\nX_{N-L} & = (f_{N-L}, \\ f_{N-L+1}, \\ f_{N-L+2}, \\ \\ldots, \\ f_{N-1} )^{\\text{T}}.\n\\end{align*}\n\nThese column vectors form the $L$-*trajectory matrix*, $\\mathbf{X}$, of the time series (hereafter just *trajectory matrix*):\n\n$$\\mathbf{X} = \\begin{bmatrix}\nf_0 & f_1 & f_2 & f_3 &\\ldots & f_{N-L} \\\\ \nf_1 & f_2 & f_3 & f_4 &\\ldots & f_{N-L+1} \\\\\nf_2 & f_3 & f_4 & f_5 &\\ldots & f_{N-L+2} \\\\\n\\vdots & \\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\nf_{L-1} & f_{L} & f_{L+1} & f_{L+2} & \\ldots & f_{N-1} \\\\ \n\\end{bmatrix}$$\n\nFrom writing out the matrix above, it is clear that the elements of the *anti-diagonals* (that is, the diagonals running from bottom-left to top-right) are equal. This type of matrix is known as a ***Hankel*** matrix.\n\nFor our toy time series, I'll set the window length to 70, and defer discussion on how to select an appropriate window length. Let $K = N - L + 1$ represent the number of columns in the trajectory matrix. **We'll refer to the columns of $\\mathbf{X}$ as the $L$-lagged vectors, and the rows as $K$-lagged vectors.**"},{"metadata":{"collapsed":true,"_cell_guid":"e0b2bbe2-3618-41c0-8c7c-8593fab49718","_uuid":"21b57609843f24bcddb0a430c852bd658f00f44d","trusted":false},"cell_type":"code","source":"L = 70 # The window length.\nK = N - L + 1 # The number of columns in the trajectory matrix.\n# Create the trajectory matrix by pulling the relevant subseries of F, and stacking them as columns.\nX = np.column_stack([F[i:i+L] for i in range(0,K)])\n# Note: the i+L above gives us up to i+L-1, as numpy array upper bounds are exclusive. ","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"b1063309-b8ad-48c4-a9ce-c68f930411e2","_uuid":"64d53a0dda38571458a44a5f6b906a1cffd1ca67"},"cell_type":"markdown","source":"Let's take a moment to appreciate the trajectory matrix in all its anti-diagonal glory:"},{"metadata":{"collapsed":true,"_cell_guid":"76351deb-32c0-44c2-852a-8e330518fdf5","_uuid":"0112dc80597ab9fc911b666274e6371694616e1b","trusted":false},"cell_type":"code","source":"ax = plt.matshow(X)\nplt.xlabel(\"$L$-Lagged Vectors\")\nplt.ylabel(\"$K$-Lagged Vectors\")\nplt.colorbar(ax.colorbar, fraction=0.025)\nax.colorbar.set_label(\"$F(t)$\")\nplt.title(\"The Trajectory Matrix for the Toy Time Series\");","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"a3f88762-ac7d-4325-87f4-ba7b2d95673a","_uuid":"a5e52e72626c19c97eba340a95c255892a762165"},"cell_type":"markdown","source":"## 2.2 Decomposing the Trajectory Matrix<a name=\"Section2.2\"></a>\nThe second step is decomposing the trajectory matrix with a [singular-value decomposition (SVD)](https://en.wikipedia.org/wiki/Singular-value_decomposition),\n\n$$\\mathbf{X} = \\mathbf{U\\Sigma V}^{\\text{T}}$$ \nwhere:\n* $\\mathbf{U}$ is an $L \\times L$ unitary matrix containing the orthonormal set of ***left singular vectors*** of $\\mathbf{X}$ as columns;\n* $\\mathbf{\\Sigma}$ is an $L \\times K$ rectangular diagonal matrix containing $L$ ***singular values*** of $\\mathbf{X}$, ordered from largest to smallest; and\n* $\\mathbf{V}$ is a $K \\times K$ unitary matrix containing the orthonormal set of ***right singular vectors*** of $\\mathbf{X}$ as columns.\n\nThe SVD of the trajectory matrix can be rewritten as \n\\begin{align*}\n    \\mathbf{X} & = \\sum_{i=0}^{d-1}\\sigma_i U_i V_i^{\\text{T}} \\\\\n               & \\equiv \\sum_{i=0}^{d-1}\\mathbf{X}_i\n\\end{align*}\nwhere $\\sigma_i$ is the $i$th singular value, $U_i$ and $V_i$ are vectors representing the $i$th columns of $\\mathbf{U}$ and $\\mathbf{V}$, respectively, $d \\le L$ is the *rank* of the trajectory matrix (which I shall elaborate on shortly), and $\\mathbf{X}_i = \\sigma_i U_i V_i^{\\text{T}}$ is the $i$th ***elementary matrix*** of $\\mathbf{X}$. The collection $\\{U_i, \\sigma_i, V_i\\}$ will be denoted the $i$th ***eigentriple*** of the SVD.\n\nTo build a picture of what all of this means, let's inspect the $\\mathbf{U}$, $\\mathbf{V}$ and $\\mathbf{\\Sigma}$ matrices in turn. \n\n### 2.2.1 The $\\mathbf{U}$ Matrix<a name=\"Section2.2.1\"></a>\n$\\mathbf{U}$ is an $L \\times L$ matrix whose columns are orthonormal, that is\n$$\n    U_i \\cdot U_j = \\left\\{\n  \\begin{array}{lr}\n    1 \\ & i = j \\\\\n    0 \\ & i \\ne j\n  \\end{array}\n\\right.\n$$\nThis means that $\\mathbf{UU}^{\\text{T}} = \\mathbf{U}^{\\text{T}}\\mathbf{U} = \\mathbf{1}$, making $\\mathbf{U}$ a unitary matrix. \n\nTo elucidate the role that $\\mathbf{U}$ plays in the expansion for $\\mathbf{X}$ above, let $Z_i = \\sigma_i V_i$ be a column vector, so that\n$$\\mathbf{X} = \\sum_{i=0}^{d-1} U_i Z_i^{\\text{T}}$$\nand each $L$-lagged column vector, $X_j$, is then given by\n$$X_j = \\sum_{i=0}^{d-1}z_{j,i}U_i$$\nwhere $z_{j,i}$ is the $j$th component of the vector $Z_i$. The expression for $X_j$ suggests that $\\mathcal{U} = \\{U_0,  \n \\ldots, U_{d-1} \\}$ is a basis set spanning the *column space* of the trajectory matrix, and $z_{j,i}$ is the $i$th coefficient of the lagged vector $X_j$ represented in the basis $\\mathcal{U}$. **In other words, the columns of the $\\mathbf{U}$ matrix form an orthonormal basis set that describes the time subseries $\\left\\{ f_i, \\ldots, f_{i+L-1}\\right\\}_{i=0}^{N-L}$ in the columns of the trajectory matrix.**  \n\n### 2.2.2 The $\\mathbf{V}$ Matrix<a name=\"Section2.2.2\"></a>\nThe matrix $\\mathbf{V}$—appearing (by convention) as its transpose in the expression for the SVD of $\\mathbf{X}$—is a $K \\times K$ matrix with orthonormal columns, which, like the $\\mathbf{U}$ matrix, makes it unitary.\n\nTo interpret the columns of $\\mathbf{V}$ in the SVD of the trajectory matrix, we first note that for any appropriately shaped matrices $\\mathbf{A}$ and $\\mathbf{B}$, $\\left(\\mathbf{AB}\\right)^{\\text{T}} = \\mathbf{B}^{\\text{T}}\\mathbf{A}^{\\text{T}}$. Taking the transpose of $\\mathbf{X}$, we therefore have\n\\begin{align*}\n\\mathbf{X}^{\\text{T}} & = \\mathbf{V \\Sigma}^{\\text{T}}\\mathbf{U}^{\\text{T}} \\\\\n                      & = \\sum_{i=0}^{d-1}V_i Y_i^{\\text{T}}\n\\end{align*}\nwhere we have set $Y_i = \\sigma_i U_i$. Then,\n$$X^{(\\text{T})}_j = \\sum_{i=0}^{d-1}y_{j,i}V_i$$ \nwhere $X^{(\\text{T})}_j$ is the $j$th column of $\\mathbf{X}^{\\text{T}}$, and $y_{j,i}$ is the $j$th component of the vector $Y_i$. This expression suggests that the $\\mathcal{V} = \\{V_0, \\ldots, V_{d-1}\\}$ is a basis set spanning the column space of $\\mathbf{X}^{\\text{T}}$, and $y_{j,i}$ is the $i$th coefficient of the lagged vector $X^{(\\text{T})}_j$ represented in the basis $\\mathcal{V}$. Equivalently, $\\mathcal{V}$ is a basis set spanning the *row space* of $\\mathbf{X}$. **That is, the columns of the $\\mathbf{V}$ matrix form an orthonormal basis set that describe the time subseries $\\{ f_i, \\ldots, f_{i+N-L}\\}_{i=0}^{L-1}$ in the rows of the trajectory matrix.**\n\n\n### 2.2.3 The $\\mathbf{\\Sigma}$ Matrix<a name=\"Section2.2.3\"></a>\nThe $\\mathbf{\\Sigma}$ matrix is an $L \\times K$ rectangular diagonal matrix containing the **singular values** of $\\mathbf{X}$. The singular values are ordered from largest to smallest, i.e. $\\sigma_0 \\ge \\sigma_1 \\ge \\ldots \\ge \\sigma_{L-1} \\ge 0$. **We can interpret $\\sigma_i$ as a scaling factor that determines the relative importance of the eigentriple $(U_i, \\sigma_i, V_i)$ in the expansion $\\mathbf{X} = \\sum_{i=0}^{d-1}\\sigma_i U_i V_i^{\\text{T}}$.**\n\nThe *Frobenius norm* of $\\mathbf{X}$, $\\lvert\\lvert \\mathbf{X} \\rvert\\rvert_{\\text{F}}$, is given by\n$$\\lvert\\lvert \\mathbf{X} \\rvert\\rvert_{\\text{F}} = \\sqrt{\\sum_{j=0}^{L-1}\\sum_{k=0}^{K-1} \\lvert x_{j,k}\\rvert^2}$$\nwhere $x_{j,k}$ denotes the element in the $j$th row and $k$th column of $\\mathbf{X}$. \n\nLet's turn our attention to the elementary matrices $\\mathbf{X}_i = \\sigma_i U_i V_i^{\\text{T}}$. Now, for an outer product such as $U_i V_i^{\\text{T}}$, we have  $\\lvert \\lvert U_i V_i^{\\text{T}} \\rvert \\rvert_{\\text{F}} = \\lvert \\lvert U_i \\rvert \\rvert_{\\text{F}} \\lvert \\lvert V_i \\rvert \\rvert_{\\text{F}}$, which is simply equal to 1 due to $U_i$ and $V_i$ being normalised. From this result, it is then clear that $\\lvert\\lvert \\mathbf{X}_i \\rvert\\rvert_{\\text{F}} = \\sigma_i$. It also turns out that \n$$\\lvert\\lvert \\mathbf{X} \\rvert\\rvert_{\\text{F}}^2 = \\sum_{i=0}^{d-1} \\sigma_i^2$$\ni.e. **the squared Frobenius norm of the trajectory matrix is equal to the sum of the squared singular values. This suggests that we can take the ratio $\\sigma_i^2 / \\lvert\\lvert \\mathbf{X} \\rvert\\rvert_{\\text{F}}^2$ as a measure of the contribution that the elementary matrix $\\mathbf{X}_i$ makes in the expansion of the trajectory matrix.**\n\nFurther, if we right-multiply the original SVD of $\\mathbf{X}$ by $\\mathbf{X}^{\\text{T}}$:\n\\begin{align*}\n    \\mathbf{XX}^{\\text{T}} & = \\mathbf{U\\Sigma V}^{\\text{T}}\\mathbf{X}^{\\text{T}} \\\\\n               & = \\mathbf{U\\Sigma V}^{\\text{T}} \\mathbf{V \\Sigma}^{\\text{T}}\\mathbf{U}^{\\text{T}} \\\\\n               & = \\mathbf{U\\Sigma} \\mathbf{\\Sigma}^{\\text{T}}\\mathbf{U}^{\\text{T}}\n\\end{align*}\nLetting the square diagonal matrix $\\mathbf{\\Sigma}^2 = \\mathbf{\\Sigma \\Sigma}^{\\text{T}}$, and multiplying on the right by $\\mathbf{U}$, gives\n$$(\\mathbf{XX}^{\\text{T}})\\mathbf{U} = \\mathbf{U}\\mathbf{\\Sigma}^2$$\nwhich—given that $\\mathbf{\\Sigma}^2$ is a diagonal matrix with elements $\\sigma_i^2$—demonstrates that the columns of $\\mathbf{U}$ are eigenvectors of the matrix $\\mathbf{XX}^{\\text{T}}$, with eigenvalues $\\{\\sigma_0^2, \\ldots , \\sigma_{L-1}^2\\}$. Following a similar argument, multiplying $\\mathbf{X}$ on the left by $\\mathbf{X}^{\\text{T}}$ shows that the columns of $\\mathbf{V}$ are eigenvectors of the matrix $\\mathbf{X}^{\\text{T}}\\mathbf{X}$, also with eigenvalues $\\{\\sigma_0^2, \\ldots , \\sigma_{L-1}^2\\}$.\n\n\n### 2.2.4 The Rank of the Trajectory Matrix<a name=\"Section2.2.4\"></a>\nWe are now equipped to talk about the *rank*, $d$, of the trajectory matrix. The columns of the trajectory matrix are a sequence of $L$-lagged vectors, which span the *trajectory space* of the time series. This space will be at most $L$-dimensional; however, if the columns in $\\mathbf{X}$ are linearly dependent, then the trajectory space will have $d < L$ dimensions. This manifests as one or more zero singular values in $\\mathbf{\\Sigma}$. **The rank of $\\mathbf{X}$ is the maximum value of $i$ such that $\\sigma_i > 0$. Put another way, $d = \\text{rank}\\{\\mathbf{X}\\}$ can be regarded as the instrinsic dimensionality of the time series' trajectory space.**\n\n(Note: for noisy real-world time series data, the trajectory space is likely to have $d = L$ dimensions.)\n\nIt is important to note that each elementary matrix $\\mathbf{X}_i$ has a rank of 1, and that the matrix \n$$\\mathbf{X}^{(r)} = \\sum_{i=0}^{r} \\mathbf{X}_i$$\nfor $r < d$, is the best rank-$r$ approximation to the trajectory matrix $\\mathbf{X}$, such that $\\lvert\\lvert \\mathbf{X} - \\mathbf{X}^{(r)}\\rvert\\rvert_F$ is minimised. **That is, we can sum the first $r$ elementary matrices to obtain an optimal, lower-dimensional approximation of $\\mathbf{X}$**.\n\n### 2.2.5 Putting it all Together<a name=\"Section2.2.5\"></a>\nLet's quickly recap everything so far: we have mapped a time series $F = \\{f_0, \\ldots, f_{N-1}\\}$ to a collection of multi-dimensional lagged vectors, $X_i = (f_i, f_{i+1}, \\ldots, f_{i+L-1})^{\\text{T}}, i = 0, \\ldots, N-L$, which together comprise the columns of the trajectory matrix $\\mathbf{X}$. We then decomposed this matrix with an SVD; in doing so, we found two orthonormal basis sets, $\\mathcal{U}$ and $\\mathcal{V}$, which span the column- and row-space, respectively, of the trajectory matrix. The SVD of $\\mathbf{X}$ can be written as\n\\begin{align*}\n    \\mathbf{X} & = \\sum_{i=0}^{d-1}\\sigma_i U_i V_i^{\\text{T}} \\\\\n               & \\equiv \\sum_{i=0}^{d-1}\\mathbf{X}_i\n\\end{align*}\nwhere $\\mathbf{X}_i$ is the $i$th elementary matrix of $\\mathbf{X}$, determined by the eigentriple  $\\{U_i, \\sigma_i, V_i\\}$. The $i$th singular value, $\\sigma_i$, determines the relative contribution of $\\mathbf{X}_i$ in the expansion of $\\mathbf{X}$ above. The integer $d \\le L$ is the intrinsic dimensionality of the time series' trajectory space, and we may choose to obtain a lower-dimensional approximation of $\\mathbf{X}$ by summing only the first $r < d$ elementary matrices.\n\nMuch of what we've covered so far is general to the SVD of *any* matrix, not just the trajectory matrix of a time series. From now on, we are going to focus on reconstructing the components of a time series from its elementary matrices. In practice, the steps involved are quite simple.\n\nEnough chat—let's decompose the trajectory matrix and form its elementary matrices:"},{"metadata":{"collapsed":true,"_cell_guid":"c9e6ec3b-0d0a-4f5d-bc38-eab28428ca53","_uuid":"8b64ba1e302404659c9a1b4d9d2ba16ac0a236de","trusted":false},"cell_type":"code","source":"d = np.linalg.matrix_rank(X) # The intrinsic dimensionality of the trajectory space.\n\n# For those interested in how to code up an SVD calculation, Numerical Recipes in Fortran 77\n# has you covered: http://www.aip.de/groups/soe/local/numres/bookfpdf/f2-6.pdf\n# Thankfully, we'll leave the actual SVD calculation to NumPy.\nU, Sigma, V = np.linalg.svd(X)\nV = V.T # Note: the SVD routine returns V^T, not V, so I'll tranpose it back here. This may seem pointless, \n# but I'll treat the Python representation of V consistently with the mathematical notation in this notebook.\n\n# Calculate the elementary matrices of X, storing them in a multidimensional NumPy array.\n# This requires calculating sigma_i * U_i * (V_i)^T for each i, or sigma_i * outer_product(U_i, V_i). \n# Note that Sigma is a 1D array of singular values, instead of the full L x K diagonal matrix.\nX_elem = np.array( [Sigma[i] * np.outer(U[:,i], V[:,i]) for i in range(0,d)] )\n\n# Quick sanity check: the sum of all elementary matrices in X_elm should be equal to X, to within a \n# *very small* tolerance:\nif not np.allclose(X, X_elem.sum(axis=0), atol=1e-10):\n    print(\"WARNING: The sum of X's elementary matrices is not equal to X!\")","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"e4e9d77e-179b-4736-8eb2-c67d4a8827b4","_uuid":"d5406cb8aff503eb96f728c52da526a75068939b"},"cell_type":"markdown","source":"Let's take a peak at the first 12 elementary matrices."},{"metadata":{"collapsed":true,"_cell_guid":"a9987d0e-3c3e-4d71-800a-f9bcbde088c0","_uuid":"1cc494833b99b4d0bd3aaa2a63be9712d34ffd02","trusted":false},"cell_type":"code","source":"n = min(12, d) # In case d is less than 12 for the toy series. Say, if we were to exclude the noise component...\nfor i in range(n):\n    plt.subplot(4,4,i+1)\n    title = \"$\\mathbf{X}_{\" + str(i) + \"}$\"\n    plot_2d(X_elem[i], title)\nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"beeba1a5-9381-49c0-aae1-fcc57ac28bbb","_uuid":"860ccd970b0e42634c378549bd2b87c31d5c7462"},"cell_type":"markdown","source":"From visual inspection of the $\\mathbf{X}_i$ above, it is obvious that the elementary matrices lack the anti-diagonal structure of the trajectory matrix. Without inspecting the $U_i$ vector associated with each $\\mathbf{X}_i$, or reconstructing a time series of each component, the appearance of the $\\mathbf{X}_i$ hints at the nature of each component, be it trend, periodicity or noise. For example, the $L$- and $K$-lagged vectors in $\\mathbf{X}_0$ and $\\mathbf{X}_1$ vary relatively slowly across the matrix, suggesting that $\\mathbf{X}_0$ and $\\mathbf{X}_1$ may be associated with the overall trend in the time series. The matrices $\\mathbf{X}_2$ to $\\mathbf{X}_5$ show large checkerboard patterns, suggesting periodicity. $\\mathbf{X}_6$ may lie somewhere between periodicity and trend. The matrices of $\\mathbf{X}_7$ onwards (and all the way up to $\\mathbf{X}_{69}$) appear to alternate quickly between a few values; these elementary matrices are likely to be associated with the noise in the original time series.\n\n(Exercise: re-run the several cells above, excluding the noise from the toy time series $F$. What happens to the number of elementary matrices? What is the intrinsic dimensionality of the noiseless series' trajectory space? For a laugh, re-run the cells with just the noise component in $F$, and note the resulting elementary matrices. Knowing that the 'process' that generated this time series is random, do you think these elementary matrices individually contain any useful information?) \n\nLet's plot the relative contributions, $\\dfrac{\\sigma_i^2}{\\sum_{k=0}^{d-1} \\sigma_k^2}$, and the cumulative contributions, $\\dfrac{\\sum_{j=0}^i \\sigma_j^2}{\\sum_{k=0}^{d-1} \\sigma_k^2}$, of the first 12 elementary matrices to the trajectory matrix of the toy time series:"},{"metadata":{"collapsed":true,"_cell_guid":"7d7e9b12-dc44-4952-a101-fc1dc3242014","_uuid":"3cb853d8c7374fed8db5c93dd24fd33263355a60","trusted":false},"cell_type":"code","source":"sigma_sumsq = (Sigma**2).sum()\nfig, ax = plt.subplots(1, 2, figsize=(14,5))\nax[0].plot(Sigma**2 / sigma_sumsq * 100, lw=2.5)\nax[0].set_xlim(0,11)\nax[0].set_title(\"Relative Contribution of $\\mathbf{X}_i$ to Trajectory Matrix\")\nax[0].set_xlabel(\"$i$\")\nax[0].set_ylabel(\"Contribution (%)\")\nax[1].plot((Sigma**2).cumsum() / sigma_sumsq * 100, lw=2.5)\nax[1].set_xlim(0,11)\nax[1].set_title(\"Cumulative Contribution of $\\mathbf{X}_i$ to Trajectory Matrix\")\nax[1].set_xlabel(\"$i$\")\nax[1].set_ylabel(\"Contribution (%)\");","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"eb4910d6-0050-46fb-bde5-e24c73c93f79","_uuid":"3d9cf20b011bf5c061f92b3dd8eb462d33c5b4fe"},"cell_type":"markdown","source":"The plots above depict the relative and cumulative contributions of the first 12 $\\mathbf{X}_i$ in the expansion $\\mathbf{X} = \\sum_{i=0}^{d-1}\\mathbf{X}_i$. The elementary matrices $\\mathbf{X}_0$ and $\\mathbf{X}_1$ contribute 52% and 22%, respectively, to the expansion of  $\\mathbf{X}$. Together, the first seven elementary matrices contribute 97%. Elementary matrices that make equal contributions to the expansion (that is, $\\sigma_i \\approx \\sigma_{i+1}$) are likely to be grouped together when reconstructing the time series, and appear as \"breaks\" in the plot of relative contributions. For example, the \"breaks\" in the plot above suggest that $\\mathbf{X}_2$ and $\\mathbf{X}_3$, and $\\mathbf{X}_4$ and $\\mathbf{X}_5$, should be grouped together. \n\nIt is important to note that the elementary matrices represent an optimal (although possibly non-unique) separation of components in trajectory space: by definition, the rows and columns of one elementary matrix are orthogonal to the rows and columns of the other elementary matrices. However, this separation may not coincide with what we would consider a useful, interpretable 'component' of the time series. In fact, there are restrictions on the types of time series components that are exactly separable under this formalism.\n\n## 2.3 Reconstructing the Time Series<a name=\"Section2.3\"></a>\nSo far, we have mapped a time series $F$ to a series of $L$-lagged vectors, forming the *trajectory matrix* of $F$. We then decomposed this matrix with a singular-value decomposition, and constructed a set of *elementary matrices* which comprise the trajectory matrix. We then gave a bit of a hand-waving explanation to classify these elementary matrices as *trend*, *periodicity* and *noise*. \n\nIn a perfect world, all the components of a time series $F = \\sum_j F^{(j)}$ would be separable, and we would have grouped the the resulting elementary matrices $\\mathbf{X}_i$ appropriately, such that\n\\begin{align*}\n\\mathbf{X} &  = \\sum_{k \\in \\mathcal{S}}\\mathbf{X}_k + \\sum_{l \\in \\mathcal{T}}\\mathbf{X}_l + \\ldots \\\\\n             &  = \\sum_j \\mathbf{X}^{(j)}\n\\end{align*}\nwhere $\\mathcal{S}$ and $\\mathcal{T}$ are disjoint (i.e. non-overlapping) sets of indices, and $\\mathbf{X}^{(j)}$ is the trajectory matrix of the time series component $F^{(j)}$. In this case, each $\\mathbf{X}^{(j)}$ would have a Hankel structure like the original trajectory matrix, and construction of each $F^{(j)}$ would be simple. However, in this imperfect real world, no component trajectory matrices will have equal values on their anti-diagonals. Therefore, we seek a process to transform an elementary matrix to a Hankel matrix, and then into a time series.\n\nTo extract a time series from the elementary matrices, we'll employ ***diagonal averaging***, which defines the values of the reconstructed time series $\\tilde{F}^{(j)}$ as averages of the corresponding anti-diagonals of the matrices $\\mathbf{X}^{(j)}$. Formally, this is represented by introducing the *Hankelisation* operator, $\\hat{\\mathcal{H}}$, that acts on the $L \\times K$ matrix $\\mathbf{X}^{(j)}$ to give a Hankel matrix $\\mathbf{\\tilde{X}}^{(j)}$; that is, \n$$\\mathbf{\\tilde{X}}^{(j)} = \\hat{\\mathcal{H}}\\mathbf{X}^{(j)}$$ \nThe element $\\tilde{x}_{m,n}$ in $\\mathbf{\\tilde{X}}^{(j)}$, for $s = m+n$, is given by\n$$\n\\tilde{x}_{m,n} = \\left\\{\n  \\begin{array}{lr}\n    \\frac{1}{s+1}\\sum_{l=0}^{s} x_{l, s-l} & \\ 0 \\le s \\le L-1 \\\\\n    \\frac{1}{L-1}\\sum_{l=0}^{L-1} x_{l, s-l} & \\ L \\le s \\le K-1 \\\\\n    \\frac{1}{K+L-s-1}\\sum_{l=s-K+1}^{L} x_{l, s-l} & \\ K \\le s \\le K+L-2 \\\\\n  \\end{array}\n\\right.\n$$\n\nAt first glance, the above looks like an impenetrable soup of matrix indices. However, all it is doing is calculating the given $\\tilde{x}_{m,n}$ by averaging the rest of the elements of the anti-diagonal wherein $\\tilde{x}_{m,n}$ belongs. The number of anti-diagonal elements to sum depends on the location of $m$ and $n$ in the matrix, and hence the index soup. **In practice, we don't need the full Hankel matrix $\\mathbf{\\tilde{X}}^{(j)}$, and can cut straight to the construction of the time series $\\tilde{F}^{(j)}$. However, I've included the definition of $\\hat{\\mathcal{H}}\\mathbf{X}^{(j)}$ above to complete the mathematical exposition of SSA.**\n\nIt is important to note that $\\hat{\\mathcal{H}}$ is a linear operator, i.e.  $\\hat{\\mathcal{H}}(\\mathbf{A} + \\mathbf{B}) = \\hat{\\mathcal{H}}\\mathbf{A} + \\hat{\\mathcal{H}}\\mathbf{B}$. Then, for a trajectory matrix $\\mathbf{X}$, \n\\begin{align*}\n\\hat{\\mathcal{H}}\\mathbf{X} & = \\hat{\\mathcal{H}} \\left( \\sum_{i=0}^{d-1} \\mathbf{X}_i \\right) \\\\\n                            &  = \\sum_{i=0}^{d-1} \\hat{\\mathcal{H}} \\mathbf{X}_i \\\\\n                            &  \\equiv \\sum_{i=0}^{d-1} \\tilde{\\mathbf{X}_i}\n\\end{align*} \nAs $\\mathbf{X}$ is already a Hankel matrix, then by definition $\\hat{\\mathcal{H}}\\mathbf{X} = \\mathbf{X}$. Therefore, the trajectory matrix can be expressed in terms of its Hankelised elementary matrices:\n$$\\mathbf{X} = \\sum_{i=0}^{d-1} \\tilde{\\mathbf{X}_i}$$\nAs a time series is uniquely determined from a Hankel matrix, the expression above also defines the time series $F$ as a sum of its components $\\tilde{F}_i$. It is up to us to group these components together, and classify them as trend, periodicity or noise, and then we're free to decide how we use them. \n\nFor fun, I've implemented a Hankelisation procedure below. However, don't get too attached to it, as it's going to be thrown away and replaced with something much more efficient and Pythonic."},{"metadata":{"collapsed":true,"_cell_guid":"9fddecaa-1352-4160-8d76-ad9c3d47740d","_uuid":"19f6bfbf357e4aba99a4ccdd06650950222868e9","trusted":false},"cell_type":"code","source":"def Hankelise(X):\n    \"\"\"\n    Hankelises the matrix X, returning H(X).\n    \"\"\"\n    L, K = X.shape\n    transpose = False\n    if L > K:\n        # The Hankelisation below only works for matrices where L < K.\n        # To Hankelise a L > K matrix, first swap L and K and tranpose X.\n        # Set flag for HX to be transposed before returning. \n        X = X.T\n        L, K = K, L\n        transpose = True\n\n    HX = np.zeros((L,K))\n    \n    # I know this isn't very efficient...\n    for m in range(L):\n        for n in range(K):\n            s = m+n\n            if 0 <= s <= L-1:\n                for l in range(0,s+1):\n                    HX[m,n] += 1/(s+1)*X[l, s-l]    \n            elif L <= s <= K-1:\n                for l in range(0,L-1):\n                    HX[m,n] += 1/(L-1)*X[l, s-l]\n            elif K <= s <= K+L-2:\n                for l in range(s-K+1,L):\n                    HX[m,n] += 1/(K+L-s-1)*X[l, s-l]\n    if transpose:\n        return HX.T\n    else:\n        return HX","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"b6d8da51-d3f8-426b-918c-3dd1da7c2e99","_uuid":"cc163faf1847f6ce525acbe0bf4242856260afae"},"cell_type":"markdown","source":"We've come this far, so we might as well Hankelise a few elementary matrices of the toy time series while we're at it:"},{"metadata":{"collapsed":true,"_cell_guid":"512f2b73-232b-4219-b127-5489875621e2","_uuid":"3e17c20b210d0460438c3d000bbbd665906fd11a","trusted":false},"cell_type":"code","source":"n = min(d, 12)\nfor j in range(0,n):\n    plt.subplot(4,4,j+1)\n    title = r\"$\\tilde{\\mathbf{X}}_{\" + str(j) + \"}$\"\n    plot_2d(Hankelise(X_elem[j]), title)\nplt.tight_layout() ","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"2f44198a-c884-460a-9703-37f606fc9f20","_uuid":"f52b7a3c4642faedd46bdf7e96379a07dc8d0e1c"},"cell_type":"markdown","source":"Inspection of the Hankelised elementary matrices of the toy time series confirms our suspicians about the elementary matrices: $\\tilde{\\mathbf{X}}_0$ and $\\tilde{\\mathbf{X}}_1$ vary slowly over the whole time series, and can be grouped together as the trend component. $\\tilde{\\mathbf{X}}_2$ and $\\tilde{\\mathbf{X}}_3$ are both periodic, with the same frequency, and can be grouped as the first periodic component. $\\tilde{\\mathbf{X}}_4$ and $\\tilde{\\mathbf{X}}_5$ are also periodic, with a different frequency to $\\tilde{\\mathbf{X}}_2$ and $\\tilde{\\mathbf{X}}_3$, and will be grouped as the second periodic component. $\\tilde{\\mathbf{X}}_6$, lacking obvious periodicity, will be grouped with the trend components. We'll lump together all components from $\\tilde{\\mathbf{X}}_7$ and beyond as noise. To summarise:\n\\begin{align*}\n\\tilde{\\mathbf{X}}^{\\text{(trend)}} & = \\tilde{\\mathbf{X}}_0 + \\tilde{\\mathbf{X}}_1 + \\tilde{\\mathbf{X}}_6 \n    & \\implies &  \\tilde{F}^{\\text{(trend)}} = \\tilde{F}_0 + \\tilde{F}_1 + \\tilde{F}_6 \\\\\n\\tilde{\\mathbf{X}}^{\\text{(periodic 1)}} & = \\tilde{\\mathbf{X}}_2 + \\tilde{\\mathbf{X}}_3 \n    & \\implies & \\tilde{F}^{\\text{(periodic 1)}} = \\tilde{F}_2 + \\tilde{F}_3  \\\\\n\\tilde{\\mathbf{X}}^{\\text{(periodic 2)}} & = \\tilde{\\mathbf{X}}_4 + \\tilde{\\mathbf{X}}_5 \n    & \\implies & \\tilde{F}^{\\text{(periodic 2)}} = \\tilde{F}_4 + \\tilde{F}_5\\\\\n\\tilde{\\mathbf{X}}^{\\text{(noise)}} & = \\tilde{\\mathbf{X}}_7 + \\ldots + \\tilde{\\mathbf{X}}_{69}\n    & \\implies & \\tilde{F}^{\\text{(noise)}} = \\tilde{F}_7 + \\ldots + \\tilde{F}_{69}\n\\end{align*}\nWhile we have defined the time series component grouping in terms of Hankelised elementary matrices, we will no longer calculate the full Hankel matrix $\\tilde{\\mathbf{X}}_i$, and instead calculate $\\tilde{F}_i$ directly from $\\mathbf{X}_i$. Further, we'll do it with two lines of Python!"},{"metadata":{"collapsed":true,"_cell_guid":"8b007026-50e0-4a3e-9258-10a1056fae7c","_uuid":"2fd8f4e741ec9a9cc937123bf98ab4e5edb85e6a","trusted":false},"cell_type":"code","source":"def X_to_TS(X_i):\n    \"\"\"Averages the anti-diagonals of the given elementary matrix, X_i, and returns a time series.\"\"\"\n    # Reverse the column ordering of X_i\n    X_rev = X_i[::-1]\n    # Full credit to Mark Tolonen at https://stackoverflow.com/a/6313414 for this one:\n    return np.array([X_rev.diagonal(i).mean() for i in range(-X_i.shape[0]+1, X_i.shape[1])])","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"ba809172-2ad3-4722-b25c-b87bd8fde4d4","_uuid":"68cecf1a5c247475ccc6c1ef1954c62268f9f634"},"cell_type":"markdown","source":"Let's go ahead and construct the first 12 elementary components, $\\tilde{F}_i$, for the toy time series."},{"metadata":{"collapsed":true,"_cell_guid":"423c15fe-7120-4754-9792-96557e31aa0f","_uuid":"23c15b51d99bd97d8aca1769b3f49f361d4f299b","trusted":false},"cell_type":"code","source":"n = min(12,d) # In case of noiseless time series with d < 12.\n\n# Fiddle with colour cycle - need more colours!\nfig = plt.subplot()\ncolor_cycle = cycler(color=plt.get_cmap('tab20').colors)\nfig.axes.set_prop_cycle(color_cycle)\n\n# Convert elementary matrices straight to a time series - no need to construct any Hankel matrices.\nfor i in range(n):\n    F_i = X_to_TS(X_elem[i])\n    fig.axes.plot(t, F_i, lw=2)\n\nfig.axes.plot(t, F, alpha=1, lw=1)\nfig.set_xlabel(\"$t$\")\nfig.set_ylabel(r\"$\\tilde{F}_i(t)$\")\nlegend = [r\"$\\tilde{F}_{%s}$\" %i for i in range(n)] + [\"$F$\"]\nfig.set_title(\"The First 12 Components of the Toy Time Series\")\nfig.legend(legend, loc=(1.05,0.1));","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"53fcae8e-5b11-4b54-a8d3-869b793f0112","_uuid":"a1a566ddd8082f392a78bfc554f7e428f158909d"},"cell_type":"markdown","source":"As mentioned earlier, the elementary components separated in the time series' trajectory space may not coincide with a single, interpretable component in the time series. For example, $\\tilde{F}_0$ and $\\tilde{F}_1$ both look vaguely like the trend—are they *really* separate components? Similarly, $\\tilde{F}_2$ and $\\tilde{F}_3$ are almost identical, except near the boundaries of the time series.\n\nWe'll introduce a way to quantify which $\\tilde{F_i}$ should be grouped together, but for a moment, let's follow our instincts and apply our earlier grouping for $\\tilde{F}^{\\text{(trend)}}$,  $\\tilde{F}^{\\text{(periodic 1)}}$,  $\\tilde{F}^{\\text{(periodic 2)}}$ and  $\\tilde{F}^{\\text{(noise)}}$, and see how the SSA-separated components compare with the original components that compose the toy time series:"},{"metadata":{"collapsed":true,"_cell_guid":"ce52e896-658c-4968-8a94-798a9aaae9d3","_uuid":"3bca431e9c6bee61224aacd3a01977e59798699b","trusted":false},"cell_type":"code","source":"# Assemble the grouped components of the time series.\nF_trend = X_to_TS(X_elem[[0,1,6]].sum(axis=0))\nF_periodic1 = X_to_TS(X_elem[[2,3]].sum(axis=0))\nF_periodic2 = X_to_TS(X_elem[[4,5]].sum(axis=0))\nF_noise = X_to_TS(X_elem[7:].sum(axis=0))\n\n# Plot the toy time series and its separated components on a single plot.\nplt.plot(t,F, lw=1)\nplt.plot(t, F_trend)\nplt.plot(t, F_periodic1)\nplt.plot(t, F_periodic2)\nplt.plot(t, F_noise, alpha=0.5)\nplt.xlabel(\"$t$\")\nplt.ylabel(r\"$\\tilde{F}^{(j)}$\")\ngroups = [\"trend\", \"periodic 1\", \"periodic 2\", \"noise\"]\nlegend = [\"$F$\"] + [r\"$\\tilde{F}^{(\\mathrm{%s})}$\"%group for group in groups]\nplt.legend(legend)\nplt.title(\"Grouped Time Series Components\")\nplt.show()\n\n# A list of tuples so we can create the next plot with a loop.\ncomponents = [(\"Trend\", trend, F_trend), \n              (\"Periodic 1\", periodic1, F_periodic1),\n              (\"Periodic 2\", periodic2, F_periodic2),\n              (\"Noise\", noise, F_noise)]\n\n# Plot the separated components and original components together.\nfig = plt.figure()\nn=1\nfor name, orig_comp, ssa_comp in components:\n    ax = fig.add_subplot(2,2,n)\n    ax.plot(t, orig_comp, linestyle=\"--\", lw=2.5, alpha=0.7)\n    ax.plot(t, ssa_comp)\n    ax.set_title(name, fontsize=16)\n    ax.set_xticks([])\n    n += 1\n\nfig.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"114b87c1-2bb5-4b93-b9a0-43c27c66fec2","_uuid":"540f11dc181c477f7dbc7ee8f8ac605ada196073"},"cell_type":"markdown","source":"After we've grouped the elementary components together, it looks like SSA has done a great job separating the original components of the toy time series—especially the two periodic components with differing frequencies and amplitudes. However, the separation isn't perfect: all the components deteriorate near the boundaries, especially the trend and second periodic component. This is common in SSA, and arises from the fact that, under the SSA formalism, most types of series (i.e. polynomial, sine, exponential, etc.) are not exactly separable. Therefore, our attempt to recover the *exact* parabolic trend and periodic components from the toy series was always doomed to fail. \n\nHowever, that does not mean components of a time series cannot be *approximately separable*, as we witnessed above. (There is also the concept of *asymptotic separability*—when the length of the time series approaches infinity—however for practical reasons we will not consider that case here.)\n\n# 3. Time Series Component Separation and Grouping<a name=\"Section3\"></a>\n\nSo far, we have grouped the eigentriples/components of the toy time series together by visual inspection; that is, we decided which components belonged together by their appearance. This is fine for a short and simple time series, however, for longer and more complicated time series, we seek a method that quantifies whether a reconstructed component $\\tilde{F}_i$ can be considered separate from another component $\\tilde{F}_j$, so we don't need to make grouping decisions by visually inspecting each $\\tilde{F}_i$.\n\nFor two reconstructed time series, $\\tilde{F}_i$ and $\\tilde{F}_j$, of length $N$, and a window length $L$, we define the *weighted inner product*, $(\\tilde{F}_i, \\tilde{F}_j)_w$ as:\n$$(\\tilde{F}_i, \\tilde{F}_j)_w = \\sum_{k=0}^{N-1} w_k \\tilde{f}_{i,k} \\tilde{f}_{j,k}$$\nwhere $\\tilde{f}_{i,k}$ and $\\tilde{f}_{j,k}$ are the $k$th values of $\\tilde{F}_i$ and $\\tilde{F}_j$, respectively, and $w_k$ is given by\n$$\nw_{k} = \\left\\{\n  \\begin{array}{lr}\n    k+1 & \\ 0 \\le k \\le L-1 \\\\\n    L & \\ L \\le k \\le K-1 \\\\\n    N - k & \\ K \\le k \\le N-1 \\\\\n  \\end{array}\n\\right.\n$$\nremembering that $K = N - L + 1$. The weight $w_k$ simply reflects the number of times $\\tilde{f}_{i,k}$ and $\\tilde{f}_{j,k}$ appear in the Hankelised matrices $\\mathbf{\\tilde{X}}_i$ and $\\mathbf{\\tilde{X}}_j$, from which the time series $\\tilde{F}_i$ and $\\tilde{F}_j$ have been obtained.\n\nPut simply, if $(\\tilde{F}_i, \\tilde{F}_j)_w = 0$, $\\tilde{F}_i$ and $\\tilde{F}_j$ are *w-orthogonal* and the time series components are separable. Of course, total w-orthogonality does not occur in real life, so instead we define a $d \\times d$ ***weighted correlation*** matrix, $\\mathbf{W}_{\\text{corr}}$, which measures the deviation of the components $\\tilde{F}_i$ and $\\tilde{F}_j$ from w-orthogonality. The elements of $\\mathbf{W}_{\\text{corr}}$ are given by\n$$W_{i,j} = \\frac{(\\tilde{F}_i, \\tilde{F}_j)_w}{\\lVert \\tilde{F}_i \\rVert_w \\lVert \\tilde{F}_j \\rVert_w}$$\nwhere $\\lVert \\tilde{F}_k \\rVert_w = \\sqrt{(\\tilde{F}_k, \\tilde{F}_k)_w}$ for $k = i,j$. The interpretation of $W_{i,j}$ is straightforward: if $\\tilde{F}_i$ and $\\tilde{F}_j$ are arbitrarily close together (but not identical), then $(\\tilde{F}_i, \\tilde{F}_j)_w \\rightarrow \\lVert \\tilde{F}_i \\rVert_w \\lVert \\tilde{F}_j \\rVert_w$ and therefore $W_{i,j} \\rightarrow 1$. Of course, if $\\tilde{F}_i$ and $\\tilde{F}_j$ are w-orthogonal, then $W_{i,j} = 0$. Moderate values of $W_{i,j}$ between 0 and 1, say $W_{i,j} \\ge 0.3$, indicate components that may need to be grouped together.\n\nWithout further ado, let's construct the w-correlation matrix for the toy time series:"},{"metadata":{"collapsed":true,"_cell_guid":"15cef5bd-e178-4204-b317-15d7f11b7cf7","_uuid":"acfda6d702610470e34e8df14659d2d2a6bebb17","trusted":false},"cell_type":"code","source":"# Get the weights w first, as they'll be reused a lot.\n# Note: list(np.arange(L)+1) returns the sequence 1 to L (first line in definition of w), \n# [L]*(K-L-1) repeats L K-L-1 times (second line in w definition)\n# list(np.arange(L)+1)[::-1] reverses the first list (equivalent to the third line)\n# Add all the lists together and we have our array of weights.\nw = np.array(list(np.arange(L)+1) + [L]*(K-L-1) + list(np.arange(L)+1)[::-1])\n\n# Get all the components of the toy series, store them as columns in F_elem array.\nF_elem = np.array([X_to_TS(X_elem[i]) for i in range(d)])\n\n# Calculate the individual weighted norms, ||F_i||_w, first, then take inverse square-root so we don't have to later.\nF_wnorms = np.array([w.dot(F_elem[i]**2) for i in range(d)])\nF_wnorms = F_wnorms**-0.5\n\n# Calculate the w-corr matrix. The diagonal elements are equal to 1, so we can start with an identity matrix\n# and iterate over all pairs of i's and j's (i != j), noting that Wij = Wji.\nWcorr = np.identity(d)\nfor i in range(d):\n    for j in range(i+1,d):\n        Wcorr[i,j] = abs(w.dot(F_elem[i]*F_elem[j]) * F_wnorms[i] * F_wnorms[j])\n        Wcorr[j,i] = Wcorr[i,j]","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"3bc6574f-e001-40c6-b47d-957ba28408ae","_uuid":"df09366694c831a2ca8b930c969817263814acfe"},"cell_type":"markdown","source":"Plot the w-correlation matrix."},{"metadata":{"collapsed":true,"_cell_guid":"6ba943df-4843-47a9-ae06-7566cc9401d5","_uuid":"2dcac6a3fe47c836753ec56ebcc490823f7e6291","trusted":false},"cell_type":"code","source":"ax = plt.imshow(Wcorr)\nplt.xlabel(r\"$\\tilde{F}_i$\")\nplt.ylabel(r\"$\\tilde{F}_j$\")\nplt.colorbar(ax.colorbar, fraction=0.045)\nax.colorbar.set_label(\"$W_{ij}$\")\nplt.clim(0,1)\nplt.title(\"The W-Correlation Matrix for the Toy Time Series\");","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"f7a6cb7b-0093-4aea-8895-7cb371f219fa","_uuid":"52eeab2e7a0936643556ef4423ea8ad0a8fe89a8"},"cell_type":"markdown","source":"The structure of $\\mathbf{W}_{\\text{corr}}$ shows a lot of correlation between the time series components, particularly in the range $7 \\le i,j \\le 69$. As these were the components we classified as belonging to the noise in the time series, it is no surprise that there are non-negligible correlations between all of them; this is a natural result of the noise having no underlying structural component that can be further separated.\n\nIt is important to note that $\\mathbf{W}_{\\text{corr}}$ is roughly split into two 'blocks': $0 \\le i,j \\le 6$, and $7 \\le i,j \\le 69$. This corresponds to two main groupings: a smoothed time series (i.e. the trend plus the two periodic components), and the residual noise. Zooming into the first seven components in $\\mathbf{W}_{\\text{corr}}$:"},{"metadata":{"collapsed":true,"_cell_guid":"0602ec98-e849-48ce-8e8f-108e54c410af","_uuid":"fbbd03845413d39ebf485be8d0a86710b8e7f14b","trusted":false},"cell_type":"code","source":"ax = plt.imshow(Wcorr)\nplt.xlabel(r\"$\\tilde{F}_i$\")\nplt.ylabel(r\"$\\tilde{F}_j$\")\nplt.colorbar(ax.colorbar, fraction=0.045)\nax.colorbar.set_label(\"$W_{ij}$\")\nplt.xlim(-0.5,6.5)\nplt.ylim(6.5,-0.5)\nplt.clim(0,1)\nplt.title(r\"W-Correlation for Components 0–6\");","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"e3069e54-9070-47d3-a9b6-b339bca0c869","_uuid":"bb4a620836f1c832c5932d5f65be42617891e2b0"},"cell_type":"markdown","source":"The initial appearance-based groupings we made for the first six components are supported by the corresponding w-correlation values. $\\tilde{F}_0$ and $\\tilde{F}_1$ have $W_{0,1} = 0.40$, suggesting they should be paired. $\\tilde{F}_1$ and $\\tilde{F}_6$ also have $W_{1,6} = 0.39$, suggesting $\\tilde{F}_6$ should also be grouped with $\\tilde{F}_0$ and $\\tilde{F}_1$ as a trend component. However, $\\tilde{F}_6$ also has a slight w-correlation with $\\tilde{F}_4$ and $\\tilde{F}_5$, but since $\\tilde{F}_4$ and $\\tilde{F}_5$ have no w-correlation with $\\tilde{F}_0$ and $\\tilde{F}_1$, we choose to keep $\\tilde{F}_6$ with $\\tilde{F}_0$ and $\\tilde{F}_1$.\n\nOur prior groupings of  $\\tilde{F}^{\\text{(periodic 1)}} = \\tilde{F}_2 + \\tilde{F}_3$ and $\\tilde{F}^{\\text{(periodic 2)}} = \\tilde{F}_4 + \\tilde{F}_5$ are clearly justified by the w-correlation matrix, with $W_{2,3} = 0.98$ and $W_{4,5} = 0.98$.\n\n(Something to ponder:  the elementary components $\\tilde{F}_2$ and $\\tilde{F}_3$ for the toy time series are almost identical, except near the boundaries, where they diverge—see the plot *The First 12 Components of the Toy Time Series* further above. However, the w-correlation for $\\tilde{F}_2$ and $\\tilde{F}_3$ is almost equal to 1. Why hasn't the boundary divergence significantly impacted the w-correlation here?)\n\n# 4. A Python Class for SSA<a name=\"Section4\"></a>\n\nSo far, the lines of code written to perform SSA on a time series have been spread through a number of cells, and are far outnumbered by the lines of code written to generate pretty plots. It is time to collect the SSA code into a handy class, imaginatively named *`SSA`*, which will form the basis for the rest of this notebook. Each instance of the class will contain the decomposition of a time series for some window length $L$, and provide useful methods to analyse, plot and reconstruct the time series.\n\nTo summarise the SSA algorithm:\n1. For a time series $F = (f_0, \\ f_1, \\ldots, \\ f_{N-1})$, and a window length $L$, form the trajectory matrix $\\mathbf{X}$, with columns given by the vectors $(f_i, \\ldots, f_{L+i-1})^{\\text{T}}$, $0 \\le i \\le N-L$.\n2. Decompose $\\mathbf{X}$ with the singular value decomposition, $\\mathbf{X} = \\sum_{i=0}^{d-1}\\sigma_i U_i V^{\\text{T}}_i$.\n3. Construct the $d$ elementary matrices $\\mathbf{X}_i = \\sigma_i U_i V^{\\text{T}}_i$.\n4. Diagonally average the $\\mathbf{X}_i$ to form the elementary time series components $\\tilde{F}_i$, such that $F = \\sum_{i=0}^{d-1} \\tilde{F}_i$.\n5. Calculate and store the weighted correlation matrix, $\\mathbf{W}_{\\text{corr}}$, for the $\\tilde{F}_i$.\n\nThe task of grouping and classifying the elementary components $\\tilde{F}_i$ is left to the user."},{"metadata":{"collapsed":true,"_cell_guid":"1b8644a0-d5e9-47cd-8840-d8d0423c0852","_uuid":"fd31e71187c5d55ff333238d20343c788c1ad643","trusted":false},"cell_type":"code","source":"class SSA(object):\n    \n    __supported_types = (pd.Series, np.ndarray, list)\n    \n    def __init__(self, tseries, L, save_mem=True):\n        \"\"\"\n        Decomposes the given time series with a singular-spectrum analysis. Assumes the values of the time series are\n        recorded at equal intervals.\n        \n        Parameters\n        ----------\n        tseries : The original time series, in the form of a Pandas Series, NumPy array or list. \n        L : The window length. Must be an integer 2 <= L <= N/2, where N is the length of the time series.\n        save_mem : Conserve memory by not retaining the elementary matrices. Recommended for long time series with\n            thousands of values. Defaults to True.\n        \n        Note: Even if an NumPy array or list is used for the initial time series, all time series returned will be\n        in the form of a Pandas Series or DataFrame object.\n        \"\"\"\n        \n        # Tedious type-checking for the initial time series\n        if not isinstance(tseries, self.__supported_types):\n            raise TypeError(\"Unsupported time series object. Try Pandas Series, NumPy array or list.\")\n        \n        # Checks to save us from ourselves\n        self.N = len(tseries)\n        if not 2 <= L <= self.N/2:\n            raise ValueError(\"The window length must be in the interval [2, N/2].\")\n        \n        self.L = L\n        self.orig_TS = pd.Series(tseries)\n        self.K = self.N - self.L + 1\n        \n        # Embed the time series in a trajectory matrix\n        self.X = np.array([self.orig_TS.values[i:L+i] for i in range(0, self.K)]).T\n        \n        # Decompose the trajectory matrix\n        self.U, self.Sigma, VT = np.linalg.svd(self.X)\n        self.d = np.linalg.matrix_rank(self.X)\n        \n        self.TS_comps = np.zeros((self.N, self.d))\n        \n        if not save_mem:\n            # Construct and save all the elementary matrices\n            self.X_elem = np.array([ self.Sigma[i]*np.outer(self.U[:,i], VT[i,:]) for i in range(self.d) ])\n\n            # Diagonally average the elementary matrices, store them as columns in array.           \n            for i in range(self.d):\n                X_rev = self.X_elem[i, ::-1]\n                self.TS_comps[:,i] = [X_rev.diagonal(j).mean() for j in range(-X_rev.shape[0]+1, X_rev.shape[1])]\n            \n            self.V = VT.T\n        else:\n            # Reconstruct the elementary matrices without storing them\n            for i in range(self.d):\n                X_elem = self.Sigma[i]*np.outer(self.U[:,i], VT[i,:])\n                X_rev = X_elem[::-1]\n                self.TS_comps[:,i] = [X_rev.diagonal(j).mean() for j in range(-X_rev.shape[0]+1, X_rev.shape[1])]\n            \n            self.X_elem = \"Re-run with save_mem=False to retain the elementary matrices.\"\n            \n            # The V array may also be very large under these circumstances, so we won't keep it.\n            self.V = \"Re-run with save_mem=False to retain the V matrix.\"\n        \n        # Calculate the w-correlation matrix.\n        self.calc_wcorr()\n            \n    def components_to_df(self, n=0):\n        \"\"\"\n        Returns all the time series components in a single Pandas DataFrame object.\n        \"\"\"\n        if n > 0:\n            n = min(n, self.d)\n        else:\n            n = self.d\n        \n        # Create list of columns - call them F0, F1, F2, ...\n        cols = [\"F{}\".format(i) for i in range(n)]\n        return pd.DataFrame(self.TS_comps[:, :n], columns=cols, index=self.orig_TS.index)\n            \n    \n    def reconstruct(self, indices):\n        \"\"\"\n        Reconstructs the time series from its elementary components, using the given indices. Returns a Pandas Series\n        object with the reconstructed time series.\n        \n        Parameters\n        ----------\n        indices: An integer, list of integers or slice(n,m) object, representing the elementary components to sum.\n        \"\"\"\n        if isinstance(indices, int): indices = [indices]\n        \n        ts_vals = self.TS_comps[:,indices].sum(axis=1)\n        return pd.Series(ts_vals, index=self.orig_TS.index)\n    \n    def calc_wcorr(self):\n        \"\"\"\n        Calculates the w-correlation matrix for the time series.\n        \"\"\"\n             \n        # Calculate the weights\n        w = np.array(list(np.arange(self.L)+1) + [self.L]*(self.K-self.L-1) + list(np.arange(self.L)+1)[::-1])\n        \n        def w_inner(F_i, F_j):\n            return w.dot(F_i*F_j)\n        \n        # Calculated weighted norms, ||F_i||_w, then invert.\n        F_wnorms = np.array([w_inner(self.TS_comps[:,i], self.TS_comps[:,i]) for i in range(self.d)])\n        F_wnorms = F_wnorms**-0.5\n        \n        # Calculate Wcorr.\n        self.Wcorr = np.identity(self.d)\n        for i in range(self.d):\n            for j in range(i+1,self.d):\n                self.Wcorr[i,j] = abs(w_inner(self.TS_comps[:,i], self.TS_comps[:,j]) * F_wnorms[i] * F_wnorms[j])\n                self.Wcorr[j,i] = self.Wcorr[i,j]\n    \n    def plot_wcorr(self, min=None, max=None):\n        \"\"\"\n        Plots the w-correlation matrix for the decomposed time series.\n        \"\"\"\n        if min is None:\n            min = 0\n        if max is None:\n            max = self.d\n        \n        if self.Wcorr is None:\n            self.calc_wcorr()\n        \n        ax = plt.imshow(self.Wcorr)\n        plt.xlabel(r\"$\\tilde{F}_i$\")\n        plt.ylabel(r\"$\\tilde{F}_j$\")\n        plt.colorbar(ax.colorbar, fraction=0.045)\n        ax.colorbar.set_label(\"$W_{i,j}$\")\n        plt.clim(0,1)\n        \n        # For plotting purposes:\n        if max == self.d:\n            max_rnge = self.d-1\n        else:\n            max_rnge = max\n        \n        plt.xlim(min-0.5, max_rnge+0.5)\n        plt.ylim(max_rnge+0.5, min-0.5)\n        ","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"466ba4c9-de5e-4491-aab4-503695db73fb","_uuid":"78f5ebaaac2de2e0233691558eed106849950a8f"},"cell_type":"markdown","source":"## 5. The Window Length<a name=\"Section5\"></a>\nWe have now established the machinery to easily investigate the effect of the window length parameter, $L$, on the decomposition of our toy time series.\n\n### 5.1 $L = 2$ <a name=\"Section5.1\"></a>\nA window length of 2 may seem like a useless choice, but it's a good place to start and watch the time series get decomposed into more and more components. We'll use the handy `SSA.components_to_df()` method to return a Pandas DataFrame and plot all elementary components at once."},{"metadata":{"collapsed":true,"_cell_guid":"13ae65a2-4962-4f25-bd43-ee85057f0b02","_uuid":"9f2bcc1b9d19434615ca9bc96e81537da39e94e5","trusted":false},"cell_type":"code","source":"F_ssa_L2 = SSA(F, 2)\nF_ssa_L2.components_to_df().plot()\nF_ssa_L2.orig_TS.plot(alpha=0.4)\nplt.xlabel(\"$t$\")\nplt.ylabel(r\"$\\tilde{F}_i(t)$\")\nplt.title(r\"$L=2$ for the Toy Time Series\");","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"182dc289-07eb-4ec6-b7b9-adf37d806cad","_uuid":"adb704d935b3a4944c54652bd71d49ec7db78ed0"},"cell_type":"markdown","source":"For $L=2$ we can only expect two elementary components to be returned. Even for such a small window length, the SSA algorithm has started to separate the high-frequency noise from the series, giving us a somewhat-denoised version of the original series in the component $\\tilde{F}_0$.\n\n## 5.2 $L = 5$<a name=\"Section5.2\"></a>\nLet's go up to a window length of 5, and see what happens to the elementary components:"},{"metadata":{"collapsed":true,"_cell_guid":"5eb6dc05-8386-4280-b256-954a6a00709e","_uuid":"e56d8e8aa99a5ab166874af463d2b590d08c0b3d","trusted":false},"cell_type":"code","source":"F_ssa_L5 = SSA(F, 5)\nF_ssa_L5.components_to_df().plot()\nF_ssa_L5.orig_TS.plot(alpha=0.4)\nplt.xlabel(\"$t$\")\nplt.ylabel(r\"$\\tilde{F}_i(t)$\")\nplt.title(r\"$L=5$ for the Toy Time Series\");\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"ca4e86ea-4233-4828-8f08-913a462b9a9e","_uuid":"6521e218cc377809419a4c8cdea152254bf71017"},"cell_type":"markdown","source":"We see that $\\tilde{F}_0$ is now a well-and-truly denoised version of the original series. $\\tilde{F}_1$ is a poorly resolved periodic component, while $\\tilde{F}_2$ to $\\tilde{F}_4$ are just noise.\n\n## 5.3 $L = 20$<a name=\"Section5.3\"></a>\n\nLet's quadruple the window length, and instead of inspecting elementary components, we'll look at the resulting w-correlation matrix and make some grouping decisions first."},{"metadata":{"collapsed":true,"_cell_guid":"04559408-ffc7-46d5-a215-757a2c3ebf9b","_uuid":"b9d176caa6b67086bfa7e6cfcff5344a9cc6f515","trusted":false},"cell_type":"code","source":"F_ssa_L20 = SSA(F, 20)\nF_ssa_L20.plot_wcorr()\nplt.title(\"W-Correlation for Toy Time Series, $L=20$\");","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"ab1f13e2-4fc8-4344-a9b2-371413503db3","_uuid":"74ebcd13a317deb527969f76f5e2ed03765915d4"},"cell_type":"markdown","source":"The w-correlation matrix for $L=20$ is split (roughly) into two blocks: $\\tilde{F}_0$ to $\\tilde{F}_3$, and $\\tilde{F}_4$ to $\\tilde{F}_{19}$. Within those blocks, the size of the $W_{i,j}$ values suggest that we need to group $\\tilde{F}_1$, $\\tilde{F}_2$ and $\\tilde{F}_3$, and group all $ \\tilde{F}_4, \\ldots,\\tilde{F}_{19}$. This grouping is certainly not ideal, as $\\tilde{F}_3$ has non-negligible w-correlation with components in the second block. We'll plot our chosen component groupings, along with $\\tilde{F}_3$ on its own, and see if we're justified in our choice of grouping: "},{"metadata":{"collapsed":true,"_cell_guid":"5f118a9a-5e98-4128-a29f-630650b22ec6","_uuid":"a3d0074cd89e4f28f96cd12b1bde418ba1d353c2","trusted":false},"cell_type":"code","source":"F_ssa_L20.reconstruct(0).plot()\nF_ssa_L20.reconstruct([1,2,3]).plot()\nF_ssa_L20.reconstruct(slice(4,20)).plot()\nF_ssa_L20.reconstruct(3).plot()\nplt.xlabel(\"$t$\")\nplt.ylabel(r\"$\\tilde{F}_i(t)$\")\nplt.title(\"Component Groupings for Toy Time Series, $L=20$\");\nplt.legend([r\"$\\tilde{F}_0$\", \n            r\"$\\tilde{F}_1+\\tilde{F}_2+\\tilde{F}_3$\", \n            r\"$\\tilde{F}_4+ \\ldots + \\tilde{F}_{19}$\",\n            r\"$\\tilde{F}_3$\"]);","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"11775d31-64be-444b-9001-0cfd49754839","_uuid":"6a53898def87e4fb340413a5158b9f645bf1292a"},"cell_type":"markdown","source":"For $L = 20$ we begin to see the trend and periodic components start to take shape. The single component $\\tilde{F}_0$ looks like the parabolic trend, and the group $\\tilde{F}_1 + \\tilde{F}_2 + \\tilde{F}_3$ is a very handsome periodicity, corresponding to the sum of the two periodic components in the original definition of the toy time series. The component $\\tilde{F}_3$ is troublesome, looking like it contributes to both noise *and* periodicity. This suggests we need to increase the window length and see if we get an improved separation of noise and periodicity.\n\n## 5.4 $L = 40$<a name=\"Section5.4\"></a>\n\nOnce again, we'll double the window length and inspect the w-correlation matrix first."},{"metadata":{"collapsed":true,"_cell_guid":"96617158-4c96-47b5-a17f-66779858705f","_uuid":"5b219e14d90c0deeaafddd1c4534fdc30bb0072c","trusted":false},"cell_type":"code","source":"F_ssa_L40 = SSA(F, 40)\nF_ssa_L40.plot_wcorr()\nplt.title(\"W-Correlation for Toy Time Series, $L=40$\");","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"372e112a-1c6c-4fd7-81d6-6aeb4c2cf745","_uuid":"b3bb7a274e1b24eb6880affa7fa1b32aedbbdb68"},"cell_type":"markdown","source":"The w-correlation matrix for $L=40$ retains the two-block structure, with $\\tilde{F}_0, \\ldots, \\tilde{F}_5$ in the first block, and $\\tilde{F}_6, \\ldots, \\tilde{F}_{39}$ in the second. Let us group the components as follows:\n\\begin{align*}\n    \\tilde{F}^{(0)} & = \\tilde{F}_0 \\\\\n    \\tilde{F}^{(1)} & = \\tilde{F}_1 + \\tilde{F}_2 + \\tilde{F}_3 \\\\\n    \\tilde{F}^{(2)} & = \\tilde{F}_4 + \\tilde{F}_5 \\\\\n    \\tilde{F}^{(3)} & = \\tilde{F}_6 + \\ldots + \\tilde{F}_{39} \\\\\n\\end{align*}\nOnce again, it can be argued that this grouping is not ideal, given the non-negligible w-correlation between, for example,  $\\tilde{F}_0$ and $\\tilde{F}_1$. "},{"metadata":{"collapsed":true,"_cell_guid":"eae846f3-1a2f-4a0b-8aa0-a59354fdfac3","_uuid":"2a2243f169465ecc64b7bf4d314605ece267d791","trusted":false},"cell_type":"code","source":"F_ssa_L40.reconstruct(0).plot()\nF_ssa_L40.reconstruct([1,2,3]).plot()\nF_ssa_L40.reconstruct([4,5]).plot()\nF_ssa_L40.reconstruct(slice(6,40)).plot(alpha=0.7)\nplt.title(\"Component Groupings for Toy Time Series, $L=40$\")\nplt.xlabel(\"$t$\")\nplt.ylabel(r\"$\\tilde{F}_i(t)$\")\nplt.legend([r\"$\\tilde{{F}}^{{({0})}}$\".format(i) for i in range(4)]);","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"71a34936-b977-4e5c-b2f6-ecba07904105","_uuid":"d850019506b8f26f1f009ebdb0f125a5f00472e4"},"cell_type":"markdown","source":"Interestingly, the trend component $\\tilde{F}^{(0)}$ has started to deteriorate at $L=40$, with notable 'kinks' in the time series. The periodicity in the toy series is now separated into two periodic components of differing amplitudes and frequencies, however with significant deteriorations near the beginning and end of the time series. At this stage our window length is 20% of the length of the time series, but the poor quality of the separated components suggests that we still need to increase the window size."},{"metadata":{"_cell_guid":"7081c568-700a-4cdf-bfb7-3a28d017a178","_uuid":"b007d5b7b6b0603d2a1bea0e0d56f123f70663fc"},"cell_type":"markdown","source":"## 5.5 $L = 60$<a name=\"Section5.5\"></a>\nWe are now approaching our first choice of $L=70$, so it's worth investigating how the decomposition converges to our original results. Inspect the w-correlation matrix first:"},{"metadata":{"collapsed":true,"_cell_guid":"34b9fe45-306d-4aca-90ec-0239a872b353","_uuid":"c1f9acde1cb2f14156b3b1a1c460b78954cdfb2d","trusted":false},"cell_type":"code","source":"F_ssa_L60 = SSA(F, 60)\nF_ssa_L60.plot_wcorr()\nplt.title(\"W-Correlation for Toy Time Series, $L=60$\");","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"cbfc596a-e00a-4acc-896e-389ab36964c3","_uuid":"5d4720c2e4d7236e1a8938cac5ca89c38b0f59b5"},"cell_type":"markdown","source":"As with the original $L=70$ result, the w-correlation matrix is now composed of two separate blocks: $\\tilde{F}_0$ to $\\tilde{F}_6$, and $\\tilde{F}_7$ to $\\tilde{F}_{59}$. From experience now, it is clear that $\\tilde{F}^{\\text{(signal)}} = \\sum_{i=0}^6 \\tilde{F}_i$ will be the combined trend and periodic components ('signal'), and $\\tilde{F}^{\\text{(noise)}} = \\sum_{i=7}^{59} \\tilde{F}_i$ will be the noise:"},{"metadata":{"collapsed":true,"_cell_guid":"0a74ac54-0276-447c-8eef-884505263835","_uuid":"5725961545d6632b04febf60bc76a8a0f85bec19","trusted":false},"cell_type":"code","source":"F_ssa_L60.reconstruct(slice(0,7)).plot()\nF_ssa_L60.reconstruct(slice(7,60)).plot()\nplt.legend([r\"$\\tilde{F}^{\\mathrm{(signal)}}$\", r\"$\\tilde{F}^{\\mathrm{(noise)}}$\"])\nplt.title(\"Signal and Noise Components of Toy Time Series, $L = 60$\")\nplt.xlabel(r\"$t$\");","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"eb8dbb67-4b37-428b-971e-3844f03e9950","_uuid":"a399e845f1d8fd08da0f3c8ccc3c5289b15454c0"},"cell_type":"markdown","source":"Viewing the w-correlation matrix for components 0–6:"},{"metadata":{"collapsed":true,"_cell_guid":"7a5c1d3a-a0b3-4390-854e-f355421af55d","_uuid":"ff22e671ab302cc1858471c9afbfdda0bc182113","trusted":false},"cell_type":"code","source":"F_ssa_L60.plot_wcorr(max=6)\nplt.title(\"W-Correlation for Toy Time Series, $L=60$ (Zoomed)\");","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"9d9d6873-5792-471f-83c6-fbb9a7d70b77","_uuid":"d822f7a4bac59b16374cab69b21a8132bc0842f8"},"cell_type":"markdown","source":"To understand why there is non-negligible w-correlation between most of the first seven components, it'll be prudent to plot all of them at once:"},{"metadata":{"collapsed":true,"_cell_guid":"a80d16f2-5c9a-45b8-bce6-5d82acb6e39e","_uuid":"2a2f84105ce59b7a19c38128ba96458d46afdb01","trusted":false},"cell_type":"code","source":"F_ssa_L60.components_to_df(n=7).plot()\nplt.title(r\"The First 7 Components of the Toy Time Series, $L=60$\")\nplt.xlabel(r\"$t$\");","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"ef8540b2-5ed1-4290-84d9-c7eddf7fc89c","_uuid":"9f9bc94b4f05dde8b3c1c98a0d7b7207e6669c01"},"cell_type":"markdown","source":"It is clear from the plot above that the trend component $\\tilde{F}_1$ contains oscillations that are the same frequency as the periodic components $\\tilde{F}_2$ and $\\tilde{F}_3$, resulting in substantial w-correlations. As we've already witnessed (see [Section 2.3](#Section2.3)), increasing the window length to $L = 70$ is enough to remove these oscillations almost entirely. It is certainly possible to increase the window length further—up to $L = N/2 = 100$ in this case—however, beyond $L=70$, the time series components generated by the leading seven eigentriples are (relatively) insensitive to the window length. \n\nIn SSA, there are no hard and fast rules to setting the perfect window length, beyond $2 \\le L \\le N/2$. However, longer window lengths (up to 30–45% of the length of the time series) are sometimes required to adequately separate underlying periodicities from the overall trend. Some trial and error is needed, but it is often easy to start at a 'large-enough' window length, and work from there.\n\n(Exercise: why is the window length $L$ restricted to the range $2 \\le L \\le N/2$? Hint: consider the expression for the SVD of the transposed trajectory matrix, $\\mathbf{X}^{\\text{T}} = \\mathbf{V \\Sigma}^{\\text{T}}\\mathbf{U}^{\\text{T}}$, for a trajectory matrix when $N/2 < L < N$.)    \n"},{"metadata":{"_cell_guid":"27325c8a-27ba-4e42-a36f-87c55eb31d12","_uuid":"c8ea7d2c9d4ca8bc4b115f97c45c2ad8c751b262"},"cell_type":"markdown","source":"## 6. MotionSense: Applying SSA to Accelerometer Data<a name=\"Section6\"></a>\n\nIn this section, we'll leave the toy time series behind and apply SSA to real-world data. To this end, I have chosen the accelerometer time series recordings in the [MotionSense](https://www.kaggle.com/malekzadeh/motionsense-dataset) dataset. More specifically, I will use accelerometer readings taken from the smart phone of a study participant walking naturally. This data is suitable for the SSA treatment because walking is a rhythmic process, with periodicities that should present themselves on accelerometer recordings.\n\n### 6.1 Loading the Data<a name=\"Section6.1\"></a>\nThe time series we'll use are located in the directory `../input/A_DeviceMotion_data/A_DeviceMotion_data/wlk_15/`. The walking time series for each subject are contained in the `sub_*.csv` files. We'll load just Subject 1 for now:"},{"metadata":{"collapsed":true,"_cell_guid":"9daafdcc-eeed-4a78-98bf-5dcceb633c61","_uuid":"db7851fac811953715d6eac11a905db2ba1a2fa4","trusted":false},"cell_type":"code","source":"walk_1 = pd.read_csv(\"../input/A_DeviceMotion_data/A_DeviceMotion_data/wlk_15/sub_1.csv\")","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"cd743641-4332-4206-a6b3-bcc4dd07d419","_uuid":"b55ef99b401443ced9cc5954521f11e83edefdbc"},"cell_type":"markdown","source":"The `userAcceleration.*` fields give the components of the instantaneous acceleration vector with respect to the fixed $x$, $y$ and $z$ coordinate axes, recorded by the participant's smart phone. For our purposes, we'll calculate the Euclidean norm of this vector—that is, we'll use the magnitude of the acceleration for our time series. The sampling rate for the measurements is 50 Hz, so we'll also transform the index to units of seconds:"},{"metadata":{"collapsed":true,"_cell_guid":"b04f453e-b652-4964-993a-51720aa985ba","_uuid":"e517d21d1bcfa34a0bc42cbbe22a2e6fe4921d97","trusted":false},"cell_type":"code","source":"cols = [\"userAcceleration.x\", \"userAcceleration.y\", \"userAcceleration.z\"]\naccel_1 = (walk_1[cols]**2).sum(axis=1)**0.5\naccel_1.index /= 50","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"c5a0f922-c6db-4fea-8495-d5e3a36bc1ad","_uuid":"70936fb3d7366c7ba66d6441b614568f2bd0958a"},"cell_type":"markdown","source":"### 6.2 Plotting the Data<a name=\"Section6.2\"></a>\nPlot the time series for the magnitude of the acceleration:"},{"metadata":{"_kg_hide-output":false,"collapsed":true,"_kg_hide-input":false,"_cell_guid":"6335adc6-b174-4056-b1a5-f34f29ecb975","_uuid":"a9d15b25f8fd542e51282647a8ad209c3ad02b7c","trusted":false},"cell_type":"code","source":"accel_1.plot()\nplt.xlabel(r\"$t$ (s)\")\nplt.ylabel(\"Acceleration (G)\")\nplt.title(\"Magnitude of Acceleration Vector While Walking\");","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"2b0dc4cb-aba9-4950-bdaf-1a68056f7479","_uuid":"157994f328eb4cfa247ed6ca869f5a5a515b41ff"},"cell_type":"markdown","source":"The instantaneous acceleration over time while walking is quite periodic, as anticipated. While we won't attempt to link the regular patterns with the features of walking (i.e. foot-fall, foot-lift, leg-swing, etc.), we'll plot a few one-second cycles on top of each other to demonstrate the periodicity in the series:  "},{"metadata":{"collapsed":true,"_cell_guid":"df7ff5c6-6775-48d6-96a5-8be41a68d9aa","_uuid":"aab1c97df9ff6cc9a878fb9c504047aad3f16559","trusted":false},"cell_type":"code","source":"accel_1.loc[4:5].plot()\naccel_1.shift(-56).loc[4:5].plot()\naccel_1.shift(-112).loc[4:5].plot()\naccel_1.shift(-168).loc[4:5].plot()\n\nplt.xlabel(r\"$t$ (s)\")\nplt.ylabel(\"Acceleration (G)\")\nplt.title(\"Example of Periodicity in Walking Time Series\");","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"978d7041-5ff8-4062-b222-cbfdf52b9fa4","_uuid":"d9d85101fba18078301c6bda3922751c6177af88"},"cell_type":"markdown","source":"### 6.3 Decomposing the Time Series With SSA<a name=\"Section6.3\"></a>\nAs a demonstration, we'll pick a 15-second (750 samples) subseries from the acceleration time series, set a window length of 7 seconds (350 samples), and apply SSA: "},{"metadata":{"collapsed":true,"_cell_guid":"b9639618-0603-40d4-8c01-deab09833d74","_uuid":"a6cee1705cae965b51933318058bcc4797006a35","trusted":false},"cell_type":"code","source":"start = 10 # seconds, not samples\nend = 25 # seconds\nwindow = 350 # samples\naccel_ssa = SSA(accel_1.loc[start:end], window)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"611e721f-0559-40a8-a36c-abf0f2ebd764","_uuid":"8b825e343a489b58853c5e4ca3207964b91c3d09"},"cell_type":"markdown","source":"As always, plot the w-correlation matrix first to help orient ourselves with the separated components."},{"metadata":{"collapsed":true,"_cell_guid":"05d5b12c-3404-404b-8ccf-b4dbfd858644","_uuid":"e832333c36db8ad472f81bf97be5b4464f94ab3a","trusted":false},"cell_type":"code","source":"accel_ssa.plot_wcorr()\nplt.title(\"W-Correlation for Walking Time Series\");","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"948c6cc1-e9e0-4513-9e31-e9fce8ba1139","_uuid":"0a6ddb23b18d4541b03c2385c135610d013c2fcd"},"cell_type":"markdown","source":"Of course, with a larger window length (and therefore a large number of elementary components), such a view of the w-correlation matrix is not the most helpful. Zoom into the w-correlation matrix for the first 50 components:  "},{"metadata":{"collapsed":true,"_cell_guid":"ce83d852-2c6a-4b1d-9bfc-97f292099d1b","_uuid":"ca644279db140d40c44f4ca27551b068ad9db824","trusted":false},"cell_type":"code","source":"accel_ssa.plot_wcorr(max=49)\nplt.title(\"W-Correlation for Walking Time Series (Zoomed)\");","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"412a3d00-8d82-4afb-9338-c13f7ab3d2e0","_uuid":"bcd4e421eb98aac9c9210355b0e0d097945bb162"},"cell_type":"markdown","source":"Beginning with the first five elementary components, we'll form the following groups and plot them:\n\n\\begin{align*}\n    \\tilde{F}^{(0)} & = \\tilde{F}_0 \\\\\n    \\tilde{F}^{(1)} & = \\tilde{F}_1 + \\tilde{F}_2 \\\\\n    \\tilde{F}^{(2)} & = \\tilde{F}_3 + \\tilde{F}_4 \\\\\n\\end{align*}"},{"metadata":{"collapsed":true,"_cell_guid":"0d6ba875-0b06-4375-adfa-2e17a55d7584","_uuid":"6a71903fba7b4b7de4d5cb77772dafa877457351","trusted":false},"cell_type":"code","source":"accel_ssa.reconstruct(0).plot()\naccel_ssa.reconstruct([1,2]).plot()\naccel_ssa.reconstruct([3,4]).plot()\naccel_ssa.orig_TS.plot(alpha=0.4)\nplt.title(\"Walking Time Series: First Three Groups\")\nplt.xlabel(r\"$t$ (s)\")\nplt.ylabel(\"Acceleration (G)\")\nlegend = [r\"$\\tilde{{F}}^{{({0})}}$\".format(i) for i in range(3)] + [\"Original TS\"]\nplt.legend(legend);","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"e2feb948-2333-4ddc-a6f7-43dfa625d2e9","_uuid":"be603f8711c2506df82efafed102fbce9a8d793c"},"cell_type":"markdown","source":"And zoom in on a four-second subseries: "},{"metadata":{"collapsed":true,"_cell_guid":"e92bb78e-e3e3-47c6-ad41-02b66b767f0f","_uuid":"3bf1736e51e82b2a8425f44452f5699f15cfdb1d","trusted":false},"cell_type":"code","source":"accel_ssa.reconstruct(0).plot()\naccel_ssa.reconstruct([1,2]).plot()\naccel_ssa.reconstruct([3,4]).plot()\naccel_ssa.orig_TS.plot(alpha=0.4)\nplt.title(\"Walking Time Series: First Three Groups (Zoomed)\")\nplt.xlabel(r\"$t$ (s)\")\nplt.ylabel(\"Acceleration (G)\")\nplt.xlim(16,20)\nplt.ylim(-0.5, 2.5)\nlegend = [r\"$\\tilde{{F}}^{{({0})}}$\".format(i) for i in range(3)] + [\"Original TS\"]\nplt.legend(legend);","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"faffb478-603c-45c2-a550-118fa50409a3","_uuid":"c11bb531c63d6b4d7e6d340a8003544de03ede3c"},"cell_type":"markdown","source":"It is clear from the two plots that the component $\\tilde{F}^{(0)}$ is the trend, while $\\tilde{F}^{(1)}$ and $\\tilde{F}^{(2)}$ are two dominating periodicities. As the long-term trend is essentially flat, the subsequent components consist of deviations above and below zero. (Note: if you subtract the mean from the time series before decomposing with SSA, the trend component $\\tilde{F}_0$ disappears.)\n\nIf we decide to add together the first five elementary components, we get a picture of the smooth, low-frequency periodicity present in the series:"},{"metadata":{"collapsed":true,"_cell_guid":"507b38c6-568d-4665-979d-387edb19f039","_uuid":"8663e3e9962d41b52fd55b83735c82d74bf3a2aa","trusted":false},"cell_type":"code","source":"accel_ssa.reconstruct(slice(0,5)).plot()\naccel_ssa.orig_TS.plot(alpha=0.4)\nplt.title(\"Walking Time Series: Low-Frequency Periodicity\")\nplt.xlabel(r\"$t$ (s)\")\nplt.ylabel(\"Acceleration (G)\")\nplt.xlim(16,20)\nplt.ylim(0, 2.5);","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"2f843772-8244-44e7-9f92-d614a034be5d","_uuid":"7ed9583bcf867bfd70622f769246a5f55505bd09"},"cell_type":"markdown","source":"Closer inspection of the original time series indicates that there are higher frequency periodicities present. Let's turn our attention to the elementary components $\\tilde{F}_5$ to $\\tilde{F}_{10}$. While these components all have some w-correlation with higher $i$ components, we'll make the following grouping:\n\n\\begin{align*}\n    \\tilde{F}^{(3)} & = \\tilde{F}_5 + \\tilde{F}_6 \\\\\n    \\tilde{F}^{(4)} & = \\tilde{F}_7 + \\tilde{F}_8 \\\\\n    \\tilde{F}^{(5)} & = \\tilde{F}_9 + \\tilde{F}_{10} \\\\\n\\end{align*}\n"},{"metadata":{"collapsed":true,"_cell_guid":"04a28eff-7d8f-48b5-b642-a0ac6f5be541","_uuid":"bc63043947cfed506d0af5084697ad82421f8f5e","trusted":false},"cell_type":"code","source":"accel_ssa.reconstruct([5,6]).plot()\n(accel_ssa.reconstruct([7,8]) + 0.25).plot() # Note: adding 0.25 to the series is for visualisation purposes only!\n(accel_ssa.reconstruct([9,10]) + 0.5).plot() # Note: adding 0.5 to the series is for visualisation purposes only!\naccel_ssa.orig_TS.plot(alpha=0.4)\nplt.title(\"Walking Time Series: Next Three Groups (Zoomed)\")\nplt.xlabel(r\"$t$ (s)\")\nplt.ylabel(\"Acceleration (G)\")\nplt.xlim(16,20)\nplt.ylim(-0.2, 2.5)\nlegend = [r\"$\\tilde{{F}}^{{({0})}}$\".format(i) for i in range(3,6)] + [\"Original TS\"]\nplt.legend(legend);","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"1c728132-de44-431d-99a3-746f37430a5d","_uuid":"bf8ae317d80cfd26b7a7a66bf7d7cf400abb42a5"},"cell_type":"markdown","source":"The grouped components $\\tilde{F}_3$, $\\tilde{F}_4$ and $\\tilde{F}_5$ are three further periodicities, with differing frequencies but approximately the same amplitudes.  In time series data with a lot of high-frequency components, SSA will typically generate a number of stable harmonic components that are well-separated from each other (i.e. low w-correlation). However, without some domain knowledge of the process generating the time series itself, it is difficult to say whether these components correspond to interpretable processes that operate independently, or if the set of components should be summed and treated as a single component. SSA is blind to reality. In this case, we'll sum these components together and inspect the periodicity they represent:"},{"metadata":{"collapsed":true,"_cell_guid":"c02f32c4-e3e4-4c58-81dc-0eca851e3cfd","_uuid":"10212b8afa77f7ab914f2c8161e724555d3e943e","trusted":false},"cell_type":"code","source":"accel_ssa.reconstruct([0] + [i for i in range(5,11)]).plot()\naccel_ssa.orig_TS.plot(alpha=0.4)\nplt.title(\"Walking Time Series: High-Frequency Periodicity (Zoomed)\")\nplt.xlabel(r\"$t$ (s)\")\nplt.ylabel(\"Acceleration (G)\")\nplt.ylim(0, 2.5)\nplt.xlim(16,20);","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"385c557c-f457-4399-bfda-f4e46372e718","_uuid":"dc0dd5e1cceaa23ef1a327fdb43fbff21acc9bb2"},"cell_type":"markdown","source":"Combining the first 11 elementary components together, and plotting the remaining 339:"},{"metadata":{"collapsed":true,"_cell_guid":"ba586ee1-dce9-44fd-a324-f93f02b8d848","_uuid":"e9c4f7acb12a8ed4268dc7be03655f0b0ea0535b","trusted":false},"cell_type":"code","source":"accel_ssa.reconstruct(slice(0,11)).plot()\naccel_ssa.reconstruct(slice(11,350)).plot()\naccel_ssa.orig_TS.plot(alpha=0.4)\nplt.ylim(-1, 2.5)\nplt.xlim(14,20)\nplt.legend([\"First 11 Components\", \"Remaining 339 Components\", \"Original TS\"]);","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"3f37804a-c4da-4dd3-88f9-9ab974c03266","_uuid":"8ac7e06013ae1ec4e9bf3a77b81369af556f9fd9"},"cell_type":"markdown","source":"The plot above demonstrates that the sum of the first 11 components adequately capture the basic, underlying periodicity of the accelerometer time series. However, the sum of the remaining 339 components (equivalent to the residual series $F_{\\text{orig}} - \\sum_{i=0}^{10} \\tilde{F}_i$) still contains large, regular spikes spaced approximately 1 second apart, likely associated with the footfall events of walking. However, the basic form of SSA presented here does not handle these types of sharp periodicities very well, instead spreading the 'signal' from these periodicities across many elementary components. Put another way, the regular part of signal cannot be separated from the noise.  I cannot offer a solution here; instead, the problem at hand (whether it's exploratory analysis, signal extraction, data preprocessing...) will dictate how situations like this are treated."},{"metadata":{"_cell_guid":"296ca82a-2748-419f-a680-7732e46b48ff","_uuid":"4b112718832b99f5a15dda6c552e0b9af316a448"},"cell_type":"markdown","source":"### 6.4 Using SSA to Extract an Individual's Walking 'Signature' <a name=\"Section6.4\"></a>\nPretend we were given a dataset consisting of unlabelled walking accelerometer time series—some of them from the same person, perhaps—and tasked with determining the number of individuals in the dataset. Let us suppose that everyone's walk has a unique 'signature' that can be extracted from accelerometer data. This signature could be compared across time series to determine whether the time series originates from the same person. Given the noise and complexity of the accelerometer time series, we therefore seek a way of extracting a simple, robust signature from the time series data.\n\nThe accelerometer time series we've seen so far contains a simple, regular low frequency component that can easily be extracted by SSA. The walking time series in the MotionSense dataset includes 24 different participants, walking on two separate occasions. We'll attempt to extract the 'signature' component from a handful of these series, and assess if it could act as an identifier for an individual. \n\n(This application of SSA has been inspired by the papers, [Persistent-homology-based gait recognition](https://arxiv.org/pdf/1707.06982.pdf) and, of course, the source of this dataset, [Protecting Sensory Data against Sensitive Inferences](https://arxiv.org/abs/1802.07802).) \n\n#### 6.4.1 Data Extraction and Processing<a name=\"Section6.4.1\"></a>\nBelow is a subroutine that will extract and return the accelerometer time series for a given subject. The parameter, `walk`  (`= 1` or `2`), specifies which of the first or second walking time series should be loaded. The time series values are converted to z-scores.\n\nSSA is then performed on the time series, using the same 15 second subseries and 7 second window length as before. The `SSA` object is then returned."},{"metadata":{"collapsed":true,"_cell_guid":"eb4de348-e79d-449a-81ed-6d33114d9699","_uuid":"7d9ed0cbfdc4abf9ddd1599a16730fe2d180156c","trusted":false},"cell_type":"code","source":"def extract_series(subject, walk):\n    walk_dict = { 1 : \"wlk_8\", 2 : \"wlk_15\"}\n    file_path = \"../input/A_DeviceMotion_data/A_DeviceMotion_data/{}/\".format(walk_dict[walk])\n    file = \"sub_{}.csv\".format(subject)\n    file_path += file\n    walk = pd.read_csv(file_path)\n    cols = [\"userAcceleration.x\", \"userAcceleration.y\", \"userAcceleration.z\"]\n    accel = (walk[cols]**2).sum(axis=1)**0.5\n    \n    # Convert acceleration values to z-scores.\n    accel = (accel - accel.mean())/accel.std()\n    \n    # Convert index to units of seconds\n    accel.index /= 50\n    \n    # Perform the SSA\n    start = 10 # seconds, not samples\n    end = 25 # seconds\n    window = 350 # samples\n    return SSA(accel.loc[start:end], window)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"1389f800-989c-44f2-9663-601bfc3562bf","_uuid":"ab4c4f4f504b9c3b0de2bb6e094f14b112596605"},"cell_type":"markdown","source":"#### 6.4.2 Comparison of Walk 'Signatures'<a name=\"Section6.4.2\"></a>\nFor this task, we'll only do an informal comparison of the walks from several participants. A more rigorous investigation would properly quantify time series similarity, ensure that the different time series are aligned as best as possible before comparison, and account for differences in walking pace. Here, we'll plot a few time series, compare them by eye, and call it a day. We will assume the low frequency components of interest are all contained in the first four eigentriples.\n\nBased on the provided participants' characteristics [here](https://github.com/mmalekzadeh/motion-sense), the participants we'll examine are as follows:\n* `3` and `5`—both females with almost-identical weights and heights;\n* `21`—male with similar age, height and weight to participants `3` and `5`; and\n* `1` and `22`—two males with similar heights and weights, approximately twice as heavy as participants `3`, `5` and `21`.\n"},{"metadata":{"_uuid":"21ac6c5a4f681278153dba2efc2161ec4b31215b","_cell_guid":"d197aa73-7704-4e07-a8e3-418f00c31e1d","scrolled":true,"trusted":false,"collapsed":true},"cell_type":"code","source":"subjects = [3,5,21,1,22]\nsubject_ssa = {}\n\n# Load and process the time series for each participant, storing the SSA objects.\n# This will take a minute or two.\nfor subject in subjects:\n    print(\"Subject #{}\".format(subject))\n    subject_ssa[subject] = extract_series(subject, 1), extract_series(subject, 2)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"3a759521-de1d-4c46-9a1c-0ce7be28ea51","_uuid":"fb3306234741f9fb980ba44209d16a0d1101f2c9"},"cell_type":"markdown","source":"Plot all of the extracted 'signatures', using the first four eigentriples of the SSA decomposition:"},{"metadata":{"_uuid":"702ca0d88ed8a182a5c58f553ca5bb4121984468","_cell_guid":"48039988-4565-4238-96e1-e199ba3889df","scrolled":false,"trusted":false,"collapsed":true},"cell_type":"code","source":"genders = [\"F\"]*2 + [\"M\"]*3\nfor subject, gender in zip(subjects, genders):\n    ssa_1, ssa_2 = subject_ssa[subject]\n    ssa_1.reconstruct(slice(0,4)).plot(color='red', figsize=(7,6))\n    ssa_2.reconstruct(slice(0,4)).plot(color='blue')\n    plt.xlim(16,20)\n    plt.ylim(-1.5,1.5)\n    plt.title(\"Participant #{} Walk 'Signature' ({})\".format(subject, gender))\n    plt.legend([\"Walk #1\", \"Walk #2\"])\n    plt.xlabel(r\"$t$ (s)\")\n    plt.ylabel(\"Z-Score\")\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"aedd340b-4797-44ab-8107-19ab43035cde","_uuid":"049e77dea350762beb44f52fb068730d7cc41606"},"cell_type":"markdown","source":"With the exception of Participant 3, the 'signature' extracted from the participants' walks are almost identical between Walk 1 and Walk 2, their relative offsets notwithstanding. While the overall time series for the male participants are very similar, we could derive some form of similarity measure to discriminate between them. (Serrà and Arcos provide a review of time series similarity measures [here](https://arxiv.org/pdf/1401.3973.pdf)).\n\nWhile this is a very rudimentary approach to discriminating between identities in time series data, it demonstrates how SSA can be easily applied to extract simple features from complicated, noisy time series.\n\n# 7. Some Final Words<a name=\"Section7\"></a>\nThrough a relatively straightforward process of embedding, decomposition and reconstruction, the technique of **singular-spectrum analysis** can extract the trend of a time series, separate underlying periodicities and remove noise. It can be used as an exploratory tool, or in the context of a more detailed analysis. \n\nThis kernel has presented both the theory and implementation of singular-spectrum analysis, and applied it to a number of time series. Please feel free to use and extend the `SSA` Python class ([Section 4](#Section4)) in other Kaggle kernels and projects—just provide a note and link to this kernel as a source.\n\nThe version of singular-spectrum analysis presented here is typically termed *basic SSA*, as a number of variants and extensions to the method have been developed. Singular-spectrum analysis can also be used for forecasting, and detecting structural changes in a time series—that is, detecting where a time series has been perturbed and taken on a new 'behaviour'. Multivariate and 2D versions of singular-spectrum analysis also exist... unfortunately, this notebook is already too long as it is.\n\nIf you have read this far, I thank you for your persistence. I hope this kernel has done justice to this wonderful method of time series analysis."}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.5","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}