This time, I further pursued the previous time series decomposition.
That is, the model is as follows
・ Decompose in R ・ Decompose in python
Almost the same as last time, you can specify seasonal data with the notation decompose_IIP $ seasonal.
data <- read.csv("./industrial/sibou_p.csv",encoding = "utf-8")
IIP <- ts(data,start=c(2008),frequency=12)
w_list <- c("seasonal","random","trend","observed")
decompose_IIP <- decompose(IIP)
for (i in w_list){
if (i=="seasonal"){
plot(decompose_IIP$seasonal)
}else if(i=="random"){
plot(decompose_IIP$random)
}else if(i=="trend"){
plot(decompose_IIP$trend)
}else if(i=="observed"){
plot(IIP)
}
}
seasonal | random |
---|---|
trend | observed |
Another code First, draw the original data. Next, since this data has a 12-month cycle, we will take a 12 moving average as a trend curve. Then, superimpose it on the above original data and draw lines.
#install.packages("forecast")
library(forecast)
plot(as.ts(IIP))
trend_beer = ma(IIP, order = 12, centre = T) #4
lines(trend_beer)
plot only draws that trend curve.
plot(as.ts(trend_beer))
Plot by subtracting the trend data from the original data.
detrend_beer = IIP - trend_beer
plot(as.ts(detrend_beer))
Then, the data with the noise minus the trend is divided into 12 pieces of data, and the average of them is calculated as the seasonal variation. 13 of them are drawn in a row. 【reference】 -Row and Column Summaries ・ Matrix
m_beer = t(matrix(data = detrend_beer, nrow = 12))
seasonal_beer = colMeans(m_beer, na.rm = T)
plot(as.ts(rep(seasonal_beer,13)))
Finally, find the random variation. The result was in agreement with what we found with the decompose () function. ** I'm looking for this calculation **.
random_beer = IIP - trend_beer - seasonal_beer
plot(as.ts(random_beer))
Finally, for various uses, output to a csv file as follows.
write.csv(as.ts(rep(seasonal_beer,1300)), file = 'decompose_seasonal1300.csv')
Similar functions are available in python. The code is below.
import pandas as pd
import statsmodels.api as sm
from statsmodels.tsa.seasonal import STL
import matplotlib.pyplot as plt
# Set figure width to 12 and height to 9
plt.rcParams['figure.figsize'] = [12, 9]
df = pd.read_csv('sibou_.csv')
series = df['Price']
print(series)
cycle, trend = sm.tsa.filters.hpfilter(series, 144)
fig, ax = plt.subplots(3,1)
ax[0].plot(series)
ax[0].set_title('Price')
ax[1].plot(trend)
ax[1].set_title('Trend')
ax[2].plot(cycle)
ax[2].set_title('Cycle')
plt.show()
Another way with python Personally, if the result of decomposition is the above method, another process is required for cycle, and this seems to be more reasonable because it can process even random parts. Also, compared to R, seasonal seems to have a larger vibration width year by year (although it is natural that R disappears because this fluctuation is averaged). In some cases, this method is the most versatile. 【reference】 -Seasonal-Trend decomposition using LOESS (STL)
stl=STL(series, period=12, robust=True)
result = stl.fit()
chart = result.plot()
chart1= result.plot(observed=False, resid=False)
chart2= result.plot(trend=False, resid=False)
plt.show()
Observed & Trend & Seasonal & Residual Trend & Seasonal Obeserved & Seasonal Also, if you do the following, you can output to one graph as a separate graph.
pl1 = result.observed
pl2 = result.trend
pl3 = result.seasonal
pl4 = result.resid
plt.plot(pl1)
plt.plot(pl2)
plt.plot(pl3)
plt.plot(pl4)
plt.show()
plt.close()
You can easily output by doing the following.
pl5 = result.seasonal.plot()
plt.show()
plt.close()
...
seasonal | resid |
---|---|
trend | observed |
And the csv file output is done as follows.
import csv
with open('sample_pl2.csv', 'w') as f:
writer = csv.writer(f, delimiter='\n')
writer.writerow(pl2)
・ I tried to arrange the method of decompose () with R and python ・ The results of each code are slightly different, but they can be almost the same, and the periodic fluctuations with noise and trends can be extracted as pure periodic fluctuations. ・ Each element decomposed component can be written to csv data and can be used secondarily. ・ Actually, this time, taking the secular change of the number of deaths as an example [preliminary figures until March](https://www.e-stat.go.jp/stat-search/files?page=1&layout=datalist&toukei=00450011&tstat=000001028897&cycle = 1 & tclass1 = 000001053058 & tclass2 = 000001053059 & result_back = 1), but no abnormal value like 2011 is seen in that range (rather, it is decreasing)
・ Next, let's convert the periodic fluctuation element into sound and listen to it. ・ Try for aperiodic fluctuation elements (aside from the test)