%matplotlib inline 

import pandas as pd
import numpy as np

from mySSA import mySSA # private lib
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 11, 4

"""
Explanation: Singular Spectrum Analysis for NH Monthly Land Temperature
Singular Spectrum Analysis (SSA) is a powerful nonparametric method of analysis and forecasting of time series. It tries to overcome the problems of finite sample length and noisiness of sampled time series not by fitting an assumed model to the available series, but by using a data-adaptive basis set, instead of the fixed sine and cosine of the BT method.
The scope of applications of SSA is very wide, from non-parametric time series decomposition and filtration to parameter estimation and forecasting. One of the differences between SSA and the methods of traditional time series analysis is the fact that SSA and SSA-related methods can be applied to quite different and not conventional for classical time series analysis problems such as exploratory analysis for data-mining and parameter estimation in signal processing, among others. For a rigorous and advanced introduction to SSA, please refer to the review paper Ghil et al. (2002) “Advanced spectral methods for climatic time series”, Rev Geophys 40(1): 3.1–3.41. 
In this notebook, we take the datasets of GISTEMP Northern Hemisphere-mean as an example (from 1880 to present). The GISS Surface Temperature Analysis (GISTEMP) is an estimate of global surface temperature change (https://data.giss.nasa.gov/gistemp/). The data is updated around the middle of every month using current data files from NOAA GHCN v3 (meteorological stations), ERSST v5 (ocean areas), and SCAR (Antarctic stations), combined as described in our December 2010 publication (Hansen et al., 2010). The data are presented as anomalies, i.e. deviations from the corresponding 1951-1980 means.
There are many implementations of SSA such as R, Mathlab and Python. The Python package of pyssa is used In this notebook.
1. Load all needed libraries
End of explanation
"""


df = pd.read_csv(r'data\NH.Ts.csv', parse_dates=True, index_col='Year', skiprows=1)
df = df['1881':'2016']
df.drop(['SON', 'MAM', 'DJF', 'JJA', 'J-D', 'D-N'], axis=1, inplace=True)
df1 = df.apply(pd.to_numeric, errors='coerce')
df1.tail()

"""
Explanation: 2. Read data
End of explanation
"""


tmp = df1.values.astype(np.float)
tmp = tmp.reshape(tmp.shape[0]*12,1)
print(tmp.shape)

ds = pd.date_range(start='1/1/1881', end='12/1/2016', freq='MS') 
ts = pd.DataFrame(tmp, columns=['TEMP'])  
ts = ts.set_index(ds)

"""
Explanation: Convert data to monthly time series
End of explanation
"""


ts.plot()

"""
Explanation: Have a quick plot
End of explanation
"""


ssa = mySSA(ts)

K = 40
suspected_seasonality = 12

ssa.embed(embedding_dimension=40, suspected_frequency=suspected_seasonality, verbose=True)

ssa.decompose(verbose=True)

"""
Explanation: 3. Carry out SSA analysis
The general procedure in SSA is as follows:
1. Embed the time series by forming a Hankel matrix of lagged window (length K) vectors.
2. Decompose the embedded time series via Singular Value Decomposition
3. Eigentripple Grouping is the process of identifying eigenvalue-eigenvector pairs as trend, seasonal and noise
4. Reconstruct the time series from the eigenvalue-eigenvector pairs identified as trend and seasonal.  This is done through a process called diagonal averaging.
Instantiate a ssa object with the time series
End of explanation
"""


ssa.view_s_contributions()

"""
Explanation: We can view the contribution of each of the signals (corresponding to each singular value) as follows:
End of explanation
"""


ssa.view_s_contributions(adjust_scale=True)

"""
Explanation: From the above we can see that most of the variation is contained in the first signal (corresponding to the first singular vector).  The scale in the image above makes it difficult to see what's going on with the remaining singular values.  
Pass the argument  adjust_scale=True to see the eigenvalue groupings more clearly.  The scale on the y-axis is meaningless and is removed in this case.
End of explanation
"""


rcParams['figure.figsize'] = 11, 2
for i in range(5):
    ssa.view_reconstruction(ssa.Xs[i], names=i, symmetric_plots=i!=0)
rcParams['figure.figsize'] = 11, 4

"""
Explanation: Now that we have the signal components (there are 15 additive components to the signal) we can loop over that range and look at each one individually.  We only look at the first 5 elements below:
The additive signal elements are stored in the object.Xs dictionary.
End of explanation
"""


streams5 = [i for i in range(5)]
reconstructed5 = ssa.view_reconstruction(*[ssa.Xs[i] for i in streams5], names=streams5, return_df=True, plot=False)

ts_copy5 = ssa.ts.copy()
ts_copy5['Reconstruction'] = reconstructed5.Reconstruction.values
ts_copy5.plot(title='Original vs. Reconstructed Time Series[First 5 RCs]');

"""
Explanation: Let's see the result if we reconstruct the signal from these 5 streams.  We add the first 5 signals and pass the result to object.view_reconstruction.
End of explanation
"""


streams15 = [i for i in range(15)]
reconstructed15 = ssa.view_reconstruction(*[ssa.Xs[i] for i in streams15], 
                                          names=streams15, return_df=True, plot=False)
ts_copy15 = ssa.ts.copy()
ts_copy15['Reconstruction'] = reconstructed15.Reconstruction.values
ts_copy15.plot(title='Original vs. Reconstructed Time Series [First 15 RCs]');

"""
Explanation: That looks very reasonable! We can get closer to the original series by adding additional signals from the remaining singular values.  Let's use the first 15 singular values:
End of explanation
"""


ssa.forecast_recurrent(steps_ahead=120, singular_values=streams15, plot=True)

"""
Explanation: 4. Forecast
To forecast using only some of the singular values, pass the list of singular value indices to be used to the forecast_recurrent method.
Take a forecast with 10 years ahead (i.e., 120 months) as an instance.
End of explanation
"""


