<jupyter_start><jupyter_text>Clustering the keywords<jupyter_code>SL = Data.iloc[[1,3,4,6,7,12,19,20,27]].sum()
UL = Data.iloc[[2,13]].sum()
RL = Data.iloc[[8,16]].sum()

ELSI = Data.iloc[[31,32,39,40,41,42,46,47,49,50,52,57]].sum()
FE = Data.iloc[[30,36,51,61]].sum()
RG = Data.iloc[[37,44,45,48,53,54,55,58,60]].sum()
RD = Data.iloc[[34,38,43,56,59]].sum()
ClustData = pd.DataFrame([SL,UL,RL,ELSI,FE,RG,RD], index=['SL','UL','RL','ELSI','FE','RG','RD']).T
stat,p = stats.normaltest(ClustData.ELSI)
print('Statistics=%.3f, p=%.3f' % (stat,p))
alpha = 0.05
if p > alpha:
    print('Pubs for ELSI look Gaussian (failed to reject H0)')
else:
    print('Pubs for ELSI do not look Gaussian (rejected H0)')
print('Kurtosis of normal distribution: {}'. format(stats.kurtosis(ClustData.ELSI) ) )
print('Skewness of normal distribution: {}'. format(stats.skew(ClustData.ELSI) ) )
plt.figure(figsize=(14,6))
plt.subplot(1,2,1)
ClustData.ELSI.hist(bins=50)
plt.title('ELSI')
plt.subplot(1,2,2)
stats.probplot(ClustData.ELSI, plot=plt);
ClustData.ELSI.describe().T
stat,p = stats.normaltest(ClustData.SL)
print('Statistics=%.3f, p=%.3f' % (stat,p))
alpha = 0.05
if p > alpha:
    print('Pubs for SL look Gaussian (failed to reject H0)')
else:
    print('Pubs for SL do not look Gaussian (rejected H0)')
print('Kurtosis of normal distribution: {}'. format(stats.kurtosis(ClustData.SL) ) )
print('Skewness of normal distribution: {}'. format(stats.skew(ClustData.SL) ) )
plt.figure(figsize=(14,6))
plt.subplot(1,2,1)
ClustData.SL.hist(bins=50)
plt.title('SL')
plt.subplot(1,2,2)
stats.probplot(ClustData.SL, plot=plt);
ClustData.SL.describe().T<jupyter_output><empty_output><jupyter_text>Correlation matrix for the clusters (Pearson and symmetric due to multiplicative commutativity)<jupyter_code>plt.figure(figsize=(12,9))
corr = ClustData.corr()
mask = np.tril(np.ones_like(corr, dtype=np.bool))
cmap = sns.cm.rocket_r
ax = sns.heatmap(corr, mask=mask, xticklabels=corr.columns.values, cmap=cmap,
            yticklabels=corr.columns.values, annot = True, annot_kws={'size' : 12} )
bottom, top = ax.get_ylim()
ax.set_ylim(bottom + 0.5, top - 0.5)
plt.tick_params(labelsize=12)
plt.title('Correlation Matrix According to Pearson', fontsize=12)

# plt.savefig('./../Manuscripts/Images/Correlation_Matrix_Clusters', dpi=150)<jupyter_output><empty_output><jupyter_text>Autocorrelation<jupyter_code># ELSI
fig, ax = plt.subplots(2, figsize=(10,6))
ax[0] = plot_acf(ClustData.ELSI, ax=ax[0])
ax[1] = plot_pacf(ClustData.ELSI, ax=ax[1])
# SL
fig, ax = plt.subplots(2, figsize=(10,6))
ax[0] = plot_acf(ClustData.SL, ax=ax[0])
ax[1] = plot_pacf(ClustData.SL, ax=ax[1])<jupyter_output><empty_output><jupyter_text>Variational Autoregression Model  Johanson's cointegration test (if there exists a linear combination of the time-series with a difference order less that that of every individual one then the test says asserts that the time-series are cointegrated) for the clustered data<jupyter_code># We only have 22 data points. Therefore we choose to keep 17 as the training set 
# and the remainder as the validation
nobs = 4
X_train, X_test = ClustData[0:-nobs], ClustData[-nobs:]

print(X_train.shape, X_test.shape)
from statsmodels.tsa.vector_ar.vecm import coint_johansen

def cointegration_test(df, alpha=0.05): 
    """Perform Johanson's Cointegration Test and Report Summary"""
    out = coint_johansen(df,-1,5)
    d = {'0.90':0, '0.95':1, '0.99':2}
    traces = out.lr1
    cvts = out.cvt[:, d[str(1-alpha)]]
    def adjust(val, length= 6): return str(val).ljust(length)

    # Summary
    print('Name   ::  Test Stat > C(95%)    =>   Signif  \n', '--'*20)
    for col, trace, cvt in zip(df.columns, traces, cvts):
        print(adjust(col), ':: ', adjust(round(trace,2), 9), ">", adjust(cvt, 8), ' =>  ' , trace > cvt)

cointegration_test(ClustData)<jupyter_output>Name   ::  Test Stat > C(95%)    =>   Signif  
 ----------------------------------------
SL     ::  168.93    > 111.7797  =>   True
UL     ::  96.27     > 83.9383   =>   True
RL     ::  52.38     > 60.0627   =>   False
ELSI   ::  26.91     > 40.1749   =>   False
FE     ::  9.92      > 24.2761   =>   False
RG     ::  1.71      > 12.3212   =>   False
RD     ::  0.12      > 4.1296    =>   False<jupyter_text>This is an iterative step. For an autoregressional model to have predictive value, one needs to satisfy stationarity. The latter can be achieved by treating the series so as to represent a non-changing behaviour on average. Possible ways of doing this is by iteratively taking the consecutive differences (differencing) of the series until stationarity (see the upcoming test) is reached, taking iteratively the logarithm of the series (if applicable) or, accordingly, some root of the data.<jupyter_code>ClustNames = ['SL', 'UL', 'RL', 'ELSI', 'FE', 'RG', 'RD']
# transform_data = pd.DataFrame(np.log(X_train.values), 
#                               index=list(range(1997,2015)), columns=ClustNames)
# TODO Bound to the augmented Dickson-Fuller test
transform_data = X_train.diff().dropna()
transform_data = transform_data.diff().dropna()
transform_data = transform_data.diff().dropna()
transform_data = transform_data.diff().dropna()
transform_data = transform_data.diff().dropna()
transform_data = transform_data.diff().dropna()
transform_data = transform_data.diff().dropna()
transform_data.head()
transform_data.describe()<jupyter_output><empty_output><jupyter_text>Stationarity check through the augmented Dickey-Fuller test<jupyter_code>class StationarityTests:
    def __init__(self, significance=.05):
        self.SignificanceLevel = significance
        self.pValue = None
        self.isStationary = None
        
    def ADF_Stationarity_Test(self, timeseries, printResults = True):
        #Dickey-Fuller test:
        adfTest = adfuller(timeseries, autolag='AIC')

        self.pValue = adfTest[1]

        if (self.pValue<self.SignificanceLevel):
            self.isStationary = True
        else:
            self.isStationary = False

        if printResults:
            dfResults = pd.Series(adfTest[0:4], index=['ADF Test Statistic','P-Value','# Lags Used','# Observations Used'])
            #Add Critical Values
            for key,value in adfTest[4].items():
                dfResults['Critical Value (%s)'%key] = value
            print('Augmented Dickey-Fuller Test Results:')
            print(dfResults)
# Check the stationarity for all the clusters with a 5% significance level
for i in ClustNames:
    sTest = StationarityTests()
    sTest.ADF_Stationarity_Test(transform_data[i])
    print("Is the time series", i, "stationary? {0}".format(sTest.isStationary))<jupyter_output>Augmented Dickey-Fuller Test Results:
ADF Test Statistic     -9.774850e+00
P-Value                 6.978921e-17
# Lags Used             2.000000e+00
# Observations Used     8.000000e+00
Critical Value (1%)    -4.665186e+00
Critical Value (5%)    -3.367187e+00
Critical Value (10%)   -2.802961e+00
dtype: float64
Is the time series SL stationary? True
Augmented Dickey-Fuller Test Results:
ADF Test Statistic     -20.472999
P-Value                  0.000000
# Lags Used              2.000000
# Observations Used      8.000000
Critical Value (1%)     -4.665186
Critical Value (5%)     -3.367187
Critical Value (10%)    -2.802961
dtype: float64
Is the time series UL stationary? True
Augmented Dickey-Fuller Test Results:
ADF Test Statistic     -1.533276e+01
P-Value                 3.925461e-28
# Lags Used             3.000000e+00
# Observations Used     7.000000e+00
Critical Value (1%)    -4.938690e+00
Critical Value (5%)    -3.477583e+00
Critical Value (10%)   -2.843868e+00
dtype: float64
Is the [...]<jupyter_text>Plots of the clusters' time-series after their transformation. The data is sparse, making the analysis somewhat speculative. One can realize though through inspection that the SL and UL clusters fluctuate beyond a mean value, making them non-stationary.<jupyter_code>fig, axes = plt.subplots(3, 3, dpi=120, figsize=(10,6))
axes[-1,-2].axis('off')
axes[-1,-1].axis('off')

for i, ax in enumerate(axes.flatten()):
    if i <=6:
        data = transform_data[transform_data.columns[i]]
        ax.plot(data, color='red', linewidth=1)
        ax.set_title(transform_data.columns[i], fontsize=8)
        ax.tick_params(labelsize=8)
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Granger's causality test for a 5% statistical significance for the stationary clusters' time-series. The resulting matrix's elements gives the p-values of the test. Anything below 5% is considered to be Granger causal.<jupyter_code># with the available data we can only look at local causal relationships, 
# not spanning more than 2 lags-years
maxlag=2
test = 'ssr_chi2test'
def grangers_causation_matrix(data, variables, test='ssr_chi2test', verbose=False):    
    """Check Granger Causality of all possible combinations of the Time series.
    The rows are response variable, columns are predictors. The values in the table 
    are the P-Values. P-Values lesser than the significance level (0.05), imply validity of 
    the Null Hypothesis, i.e. that the coefficients of the corresponding past values are
    zero, that is, that X does not cause Y can be rejected.

    data      : pandas dataframe containing the time series variables
    variables : list containing names of the time series variables.
    """
    df = pd.DataFrame(np.zeros((len(variables), len(variables))), columns=variables, index=variables)
    for c in df.columns:
        for r in df.index:
            test_result = grangercausalitytests(data[[r, c]], maxlag=maxlag, verbose=False)
            p_values = [round(test_result[i+1][0][test][1],4) for i in range(maxlag)]
            if verbose: print(f'Y = {r}, X = {c}, P Values = {p_values}')
            min_p_value = np.min(p_values)
            df.loc[r, c] = min_p_value
    df.columns = [var + '_x' for var in variables]
    df.index = [var + '_y' for var in variables]
    return df<jupyter_output><empty_output><jupyter_text>The Granger causation matrix with a 5% stat. significance (anything below 5% is supposed to be Granger causal)<jupyter_code># StatTSeries = X_train.drop(['SL','UL'], axis=1)
mask = GClusters.values >= 0.05  # keep only values smaller than 5%
GClusters = grangers_causation_matrix(transform_data, variables=transform_data.columns)
plt.figure(figsize=(12,9))
cmap = sns.cm.rocket_r
ax = sns.heatmap(GClusters, mask = mask, xticklabels=GClusters.columns.values, cmap=cmap,
            yticklabels=GClusters.index, annot = True, annot_kws={'size' : 12}, vmin=0, vmax=0.04 )
bottom, top = ax.get_ylim()
ax.set_ylim(bottom + 0.5, top - 0.5)
ax.tick_params(labelsize=12)
plt.title('Granger causality matrix of X cluster causing Y', fontsize=14)
# plt.savefig('./../Manuscripts/Images/Granger_Causation_Clusters', dpi=150)
grangercausalitytests(transform_data[['ELSI','ELSI']], 1)
import statsmodels.tsa.api as smt
from statsmodels.tsa.api import VAR
transform_data.shape
model = VAR(transform_data)
for i in [1,2,3]:
    result = model.fit(i)
    print('Lag Order =', i)
    print('AIC : ', result.aic)
    print('BIC : ', result.bic)
    print('FPE : ', result.fpe)
    print('HQIC: ', result.hqic, '\n')<jupyter_output>Lag Order = 1<jupyter_text>Discrete Fourier Transform (determining the dominant frequencies)<jupyter_code>y = {}
x = {}
for i in ClustNames:
    y[i] = np.fft.rfft(ClustData[i])
    x[i] = ClustData[i].index.shape[0]/np.arange(1,y[i].shape[0])
    y[i] = np.abs(y[i])[1:] ** 2.
fig, axes = plt.subplots(3, 3, dpi=120, figsize=(10,6))
axes[-1,-2].axis('off')
axes[-1,-1].axis('off')

for i, ax in zip(ClustNames, axes.flatten()):
    ax.bar(x[i],y[i])
    # plt.xscale('log')
    ax.set_yscale('log')
    ax.set_ylabel('Intensity', fontsize=12)
    ax.set_xlabel('Periods (in 22 Years)', fontsize=12)
    ax.tick_params(labelsize=12)
    ax.set_title(i)

plt.tight_layout()
# plt.savefig('./../Manuscripts/Images/Cluster_DFT', dpi=150)<jupyter_output><empty_output><jupyter_text>Wavelet analysis (tentative) **WARNING**: This part of the analysis uses an [external package](https://github.com/aaren/wavelets) for wavelet analysis. This package is copyrighted by <NAME> under the MIT License, granting free use of the code.We use [Morlet](https://en.wikipedia.org/wiki/Morlet_wavelet) basis function.<jupyter_code>from wavelets import WaveletAnalysis
dt = 0.1  # sample spacing

wa = WaveletAnalysis(transform_data['ELSI'].values, dt=dt)
# wa = WaveletAnalysis(ClustData['ELSI'].values, dt=dt)

# wavelet power spectrum
power = wa.wavelet_power

# scales 
scales = wa.scales

# associated time vector
t = wa.time

# reconstruction of the original data
rx = wa.reconstruction()<jupyter_output>/home/stavrogin/anaconda3/lib/python3.7/site-packages/wavelets/transform.py:104: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  wavelet_data[slices],<jupyter_text>According to wavelet theory and the documentation of the algorithm, the time-series  must be localized in the time-frequency domain and have a zero mean. So our  transformed (static) ones best serve this property. The results are obviously radically different if we use the original time-series<jupyter_code>fig, ax = plt.subplots()
T, S = np.meshgrid(t, scales)
cax = ax.contourf(T, S, power, 100, cmap='bwr')
ax.set_yscale('log')
ax.set_xlabel('Time', fontsize=12)
ax.set_ylabel('Scale', fontsize=12)
cbar = fig.colorbar(cax)
cbar.set_label('Power (Pubs No.)', fontsize=12)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13802386495925495, 'fraction_numerical': 0.04118160651920838, 'mean_word_length': 4.186792452830189, 'pattern_counts': {'":': 0, '<': 48, '<?xml version=': 0, '>': 70, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '474113', 'license': 'MIT', 'n_tokens_mistral': 5166, 'n_tokens_neox': 4731, 'n_words': 1209, 'path': 'G_Causality.ipynb', 'repo_name': 'Yperidis/TserSUITE', 'size': '13751'}