How to plot statsmodels timeseries plots side by side and customize x axis in Python

Question

I am creating these timeseries plots specifically stl decomposition and already managed to get all the plots into one. The issue I am having is having them shown side by side like the solution [here][1]. I tried the solution on the link but it did not work, instead I kept getting an empty plot on the top. I have four time series plots and managed to get them outputted on the bottom of each other however I would like to have them side by side or two side by side and the last two on the bottom side by side.

Then for the dates on the xaxis, I have already tried using ax.xaxis.set_major_formatter(DateFormatter('%b %Y')) but it is not working on the code below since the res.plot function won't allow it.

I have already searched everywhere but I can't find the solution to my issue. I would appreciate any help.

Answer 1

Here is an example using artificial data. The main idea is to group the outputs in to DataFrames and then to plot these using the pandas plot function.

Note that I had to change your code to use stl2, stl3 and stl4 when fitting.

from statsmodels.tsa.seasonal import STL
import matplotlib.pyplot as plt
import seaborn as sns
from pandas.plotting import register_matplotlib_converters 
from matplotlib.dates import DateFormatter

register_matplotlib_converters()
sns.set(style='whitegrid', palette = sns.color_palette('winter'), rc={'axes.titlesize':17,'axes.labelsize':17, 'grid.linewidth': 0.5})
plt.rc("axes.spines", top=False, bottom = False, right=False, left=False)
plt.rc('font', size=13)
plt.rc('figure',figsize=(17,12))


idx = pd.date_range("1-1-2020", periods=200, freq="M")
seas = 10*np.sin(np.arange(200) * np.pi/12)
trend = np.arange(200) / 10.0
seatr = pd.Series(trend + seas + np.random.standard_normal(200), name="Seattle", index=idx)
latr = pd.Series(trend + seas + np.random.standard_normal(200), name="LA", index=idx)
sftr = pd.Series(trend + seas + np.random.standard_normal(200), name="SF", index=idx)
phtr = pd.Series(trend + seas + np.random.standard_normal(200), name="Philly", index=idx)

stl = STL(seatr, seasonal=13)
res = stl.fit()

stl2 = STL(latr, seasonal=13)
res2 = stl2.fit()

stl3 = STL(sftr, seasonal=13)
res3 = stl3.fit()

stl4 = STL(phtr, seasonal=13)
res4 = stl4.fit()

data = pd.concat([seatr, latr, sftr, phtr], 1)
trends = pd.concat([res.trend, res2.trend, res3.trend, res4.trend], 1)
seasonals = pd.concat([res.seasonal, res2.seasonal, res3.seasonal, res4.seasonal], 1)
resids = pd.concat([res.resid, res2.resid, res3.resid, res4.resid], 1)

fig, axes = plt.subplots(4,1)
data.plot(ax=axes[0])
trends.plot(ax=axes[1])
seasonals.plot(ax=axes[2])
resids.plot(ax=axes[3])

This produces: