Una prova per rimuovere la componente stagionale da una serie tempo usando la libreria StatsModels di Python
In alto dati originali. al centro trend, seguend dati periodici stimati ed in basso calcolo dei residui
#!/usr/bin/env python
# coding: utf-8
# In[1]:
from IPython.display import display
import numpy as np
import pandas as pd
pd.set_option('display.max_rows', 15)
pd.set_option('display.max_columns', 500)
pd.set_option('display.width', 1000)
import matplotlib.pyplot as plt
from datetime import datetime
from datetime import timedelta
from pandas.plotting import register_matplotlib_converters
from mpl_toolkits.mplot3d import Axes3D
from statsmodels.tsa.stattools import acf, pacf
from statsmodels.tsa.statespace.sarimax import SARIMAX
import statsmodels.api as sm
import statsmodels.formula.api as smf
register_matplotlib_converters()
from time import time
import seaborn as sns
sns.set(style="whitegrid")
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
from sklearn.covariance import EllipticEnvelope
import warnings
warnings.filterwarnings('ignore')
RANDOM_SEED = np.random.seed(0)
# # Carica i dati
# In[2]:
def parser(s):
return datetime.strptime(s, '%Y-%m-%d %H:%M:%S')
# In[3]:
#read data
est = pd.read_csv('/home/luca/luca/liscione.csv',sep=";",parse_dates=[0], index_col=0, date_parser=parser,header=0)
print(est)
# In[4]:
est = est.asfreq(freq='1H', method='ffill')
# ### Introduce an Anomaly
# In[5]:
idx = pd.IndexSlice
# In[6]:
start_date = datetime(2023,10,5,0,0,0)
end_date = datetime(2023,12,31,23,59,0)
lim_est = est[start_date:end_date]
# In[7]:
lim_est
# In[8]:
del lim_est['Temp']
# In[9]:
plt.figure(figsize=(10,4))
plt.plot(lim_est)
plt.title('Est', fontsize=20)
plt.ylabel('Est(mm)', fontsize=16)
for year in range(start_date.year,end_date.year):
plt.axvline(pd.to_datetime(str(year)+'-01-01'), color='k', linestyle='--', alpha=0.2)
# ## Seasonal Decompose
# In[10]:
from statsmodels.tsa.seasonal import seasonal_decompose
import matplotlib.dates as mdates
# In[18]:
plt.rc('figure',figsize=(12,8))
plt.rc('font',size=15)
result = seasonal_decompose(lim_est,model='additive')
print(result.trend)
plt.savefig('trend')
fig = result.plot()
fig.savefig('trend.png') # save the figure to file
result.trend.to_csv('/home/luca/tempo.csv', date_format='%Y-%m-%d_%H',sep=";")
# In[12]:
plt.rc('figure',figsize=(12,6))
plt.rc('font',size=15)
fig, ax = plt.subplots()
x = result.resid.index
y = result.resid.values
print(x.size)
ax.plot_date(x, y, color='black',linestyle='--')
ax.set_ylim([-0.7, 0.7])
fig.autofmt_xdate()
plt.savefig('residual')
plt.show()
# In[13]:
temp = est[start_date:end_date]
del temp['Est']
print(temp.index)
temp.index.to_csv('tempo2.csv', date_format='%Y-%m-%d_%H',delimiter=";")
#np.savetxt('tempo.csv',temp.index,delimiter=";")
# In[ ]:
# In[ ]:
fig, ax1 = plt.subplots()
color = 'tab:red'
ax1.set_xlabel('time')
ax1.set_ylabel('Temp', color=color)
ax1.plot(temp, color=color)
ax1.tick_params(axis='y', labelcolor=color)
ax2 = ax1.twinx() # instantiate a second axes that shares the same x-axis
color = 'tab:blue'
ax2.set_ylabel('Est', color=color) # we already handled the x-label with ax1
ax2.plot(result.trend, color=color)
ax2.tick_params(axis='y', labelcolor=color)
plt.savefig('sovrapposto')
plt.show()
# In[ ]:
np.savetxt('trend.csv',result.trend.values,delimiter=";")
np.savetxt('temp.csv',temp.values,delimiter=";")
# In[ ]:
import matplotlib.animation as animation
fig,ax = plt.subplots()
scat = ax.scatter(result.trend.values, temp.values)
plt.gca().invert_xaxis()
ax.set(xlabel='Estensimetro (mm)')
ax.set(ylabel='Temperatura (C)')
def update(i):
scat.set_offsets((result.trend.values[i],temp.values[i]))
return scat,
ani = animation.FuncAnimation(fig,update, repeat=True, frames=len(temp)-1, interval=10)
#writer = animation.PillowWriter(fps=15)
#ani.save('animation.gif',writer=writer)
#plt.savefig('scatterplot')
plt.show()
