import pandas as pd
from statsmodels.tsa.seasonal import seasonal_decompose
import matplotlib.pyplot as plt
from pandas.plotting import lag_plot
from PythonTsa.plot_acf_pacf import acf_pacf_fig
from statsmodels.tsa.stattools import acf
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from statsmodels.stats.diagnostic import acorr_ljungbox
x = pd.read_csv('gdpquarterlychina1992.1-2017.4.csv', header = 0) 
dates = pd.date_range(start = '1992', periods = len(x), freq = 'Q')
x.index = dates
x = pd.Series(x['GDP'])
xdeca = seasonal_decompose(x, model = 'additive')
xdecm = seasonal_decompose(x, model = 'multiplicative')
xdeca.plot();  plt.show() 
# aditivo   plt.title('Additive Decomposition')

xdecm.plot();  plt.show() 
# multiplicativo  plt.title('Multiplicative Decomposition')

fig = plt.figure()
lag_plot(xdeca.resid, ax = fig.add_subplot(211))
plt.title('Additive Decomposition')

import numpy as np
logxts = np.log(x)
dlogxts = logxts.diff(1)
dlogxts = dlogxts.dropna() #delete "NaN"
dlogxts.plot(marker = 'o', markersize = 5)
plt.title('Difference of Logarithm of the GDP China');plt.show()

acf_pacf_fig(dlogxts, both = True, lag = 17)
r,q,p = acf(dlogxts,nlags = 35,qstat = True)
#r for ACF; q for Ljung-Box statistics; p for p-values
p
## array([3.39833109e-007, 1.28194353e-006, 3.83072574e-012, 2.34070680e-032,
##        9.67850548e-037, 3.70104149e-036, 1.86041495e-041, 4.46934351e-060,
##        3.47439176e-064, 1.40602250e-063, 8.85993838e-069, 7.09214985e-086,
##        1.28800665e-089, 5.30996928e-089, 4.54253985e-094, 1.43673615e-109,
##        4.84415723e-113, 2.06493251e-112, 3.18394507e-117, 3.68768712e-131,
##        2.27821600e-134, 9.91264483e-134, 4.23843931e-138, 9.39523401e-151,
##        1.11440844e-153, 4.70240196e-153, 6.77435254e-157, 1.88549878e-168,
##        3.94424975e-171, 1.66784396e-170, 6.84520898e-174, 2.06898641e-184,
##        6.58838274e-187, 2.80936462e-186, 2.78767874e-189])
q1,p1 = acorr_ljungbox(dlogxts, lags = 35, boxpierce = False)
p1
## 'lb_pvalue'
ddlogxts = dlogxts.diff(1)
ddlogxts = ddlogxts.dropna()
acf_pacf_fig(ddlogxts, both = True, lag = 17); plt.show()

r2,q2,p2 = acf(ddlogxts,nlags = 35,qstat = True)
p2
## array([7.49714347e-013, 1.00337097e-015, 8.67676846e-026, 4.48166322e-046,
##        1.61237455e-055, 3.17115154e-058, 1.08902326e-067, 1.11290994e-086,
##        1.96671898e-095, 5.77655549e-098, 7.28187337e-107, 1.83851603e-124,
##        1.77387077e-132, 7.75972346e-135, 3.54343109e-143, 2.01700636e-159,
##        9.31002746e-167, 6.79401417e-169, 1.56784081e-176, 2.27657823e-191,
##        4.73345486e-198, 5.93576596e-200, 7.64173376e-207, 2.32055469e-220,
##        1.60912563e-226, 3.16307936e-228, 2.13290288e-234, 1.09006511e-246,
##        2.49405243e-252, 7.89815787e-254, 2.29778552e-259, 1.40236100e-270,
##        9.13886515e-276, 4.23335661e-277, 4.32339487e-282])
lag_plot(xdecm.resid, ax = fig.add_subplot(212))
plt.title('Multiplicative Decomposition')
#fig.set_size_inches(18.5, 10.5)
fig.tight_layout(pad = 1.5)

df = pd.DataFrame(data = { 'resid'  : xdecm.resid})
df = df.dropna()
acf_pacf_fig(df, both = True, lag = 20)
plt.show()

xhwfit = ExponentialSmoothing(x,trend = 'add',seasonal = 'add', seasonal_periods = 4).fit()
## C:\Users\amrof\AppData\Local\Programs\Python\Python310\lib\site-packages\statsmodels\tsa\holtwinters\model.py:915: ConvergenceWarning: Optimization failed to converge. Check mle_retvals.
##   warnings.warn(
ax1 = plt.subplot(411);x.plot();
plt.setp(ax1.get_xticklabels(),visible = False); 
plt.ylabel('GDP')
ax2 = plt.subplot(412,sharex = ax1); xhwfit.level.plot();
plt.setp(ax2.get_xticklabels(),visible = False); plt.ylabel('Level')
#ax3 = plt.subplot(413,sharex = ax1);xhwfit.slope.plot();
#plt.setp(ax3.get_xticklabels(),visible = False); plt.ylabel('Trend')

ax4 = plt.subplot(414, sharex = ax1); xhwfit.season.plot();
plt.ylabel('Season') 
fig.set_size_inches(18.5, 10.5)
fig.tight_layout(pad = 1.5); plt.show()

xhwfit.resid.plot(); plt.show()

lag_plot(xhwfit.resid); plt.show()

acf_pacf_fig(xhwfit.resid, both = True, lag = 20); plt.show()

y = xhwfit.forecast(4)
y
## 2018-03-31    202227.765000
## 2018-06-30    222930.713282
## 2018-09-30    234460.516819
## 2018-12-31    254767.460457
## Freq: Q-DEC, dtype: float64
z = xhwfit.predict(start = '1992-03-31',end = '2018-12-31')
z = pd.DataFrame(z,columns = {'Predict'})
zx = z.join(x)
Predict, = plt.plot(zx['Predict'],marker = '.',label = 'Predict')
GDP, = plt.plot(zx['GDP'], linewidth = 1.0, label = 'GDP')
plt.title('Chinese Quarterly GDP and Predict')
plt.legend(handles = [Predict, GDP])
fig.tight_layout(pad = 1.5)
r,q,p = acf(xhwfit.resid,qstat = True, nlags = 20)
p
## array([3.07406843e-03, 8.45171570e-03, 7.15196695e-04, 9.81704607e-05,
##        1.39651441e-04, 1.31680677e-04, 4.77447016e-06, 8.27957022e-06,
##        1.17272256e-06, 8.70972299e-07, 1.77801388e-06, 1.26752197e-11,
##        1.35825421e-11, 3.43336853e-11, 5.76927188e-11, 1.08373071e-12,
##        2.35375418e-12, 7.82276682e-13, 4.43388691e-15, 1.05142824e-14])