Time Series Forecast Of Sunspot Counts
Daniel Kurnia
June 1, 2021
1 Introduction
sunspot image
Sunspots, as seen in the above image are dark areas on the surface of the sun. They appear dark because those areas have lower temperature than the bright areas. The sunspots are caused by concentrated magnetic field that hinders heat convection (source). Although they are relatively cooler than the surrounding bright areas, their temperatures are still high, about thousands of degrees celsius source. The sunspots can appear and disappear, their number changed throughout the years. The sunspot count has been observed to repeat roughly on 11 years, this period is called the sunspot cycle. The appearance and dissapearance of sunspots can affect the earth’s climate according to a government website. The same site cites that from 1645 - 1715 there was a period of near zero sunspot activity which coincides with “The Little Ice Age”. The sunspot activity is proposed as one of the factors that caused this phenomenon. Readers who want to know more about the little ice age can read this wikipedia article
So, knowing the sunspot cycle may help in predicting the Earth’s climate change. Now we will use the sunspot historic data to confirm or reject this trend and to see if there are additional patterns in the sun’s sunspot count.
2 The data
The monthly mean sunspot numbers data are from Kaggle. The website shows that the time range is from January 1749 to January 2021.
3 Data Import, Wrangling and EDA
sunspot <- read.csv("Sunspots.csv")
glimpse(sunspot)#> Rows: 3,265
#> Columns: 3
#> $ X <int> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11…
#> $ Date <chr> "1749-01-31", "1749-02-28", "1749-03…
#> $ Monthly.Mean.Total.Sunspot.Number <dbl> 96.7, 104.3, 116.7, 92.8, 141.7, 139…We will convert the date type from char to date in order to apply pad to it to check if there is any missing time data. Additionally we will remove the index column.
sun_clean <- sunspot %>%
select(-X) %>%
mutate(Date = ymd(Date)) %>%
thicken(interval ="month", colname = "yearmonth") %>% # to convert "daily" observations to monthly observations as there is only one day per month
select(-"Date") %>% # remove the original date column to use the yearmonth column as the new date column as basis for padding
pad(interval = "month")
anyNA(sun_clean)#> [1] FALSEThere are no missing intervals in our data so we can proceed with time series object creation. We will first use a yearly frequency.
sunspot_ts <- ts(data=sun_clean$Monthly.Mean.Total.Sunspot.Number,
start = c(1749,1),
frequency = 12)
autoplot(sunspot_ts)
The sunspot data shows no obvious increasing and decreasing trend over more than two and a half centuries but it seems that the peaks is pretty regular which shows indication of seasonality. Now let us try to decompose the time series.
autoplot(decompose(sunspot_ts))
The trend, seasonality, and remainder all still show repeated patterns which mean that there is still a seasonality we have not captured. Let’s make the frequency to be 11 years as suggested by literature.
sunspot_11yr <- ts(data=sun_clean$Monthly.Mean.Total.Sunspot.Number,
start = c(1749,1),
frequency = 12*11)
autoplot(decompose(sunspot_11yr))
Now the seasonality is separated into definite peaks and the trend is less bumpy but it still show a kind of pattern. This confirms the observations by astronomers that the sunspot count has a seasonality that repeats roughly every 11 years. The remainder also seems to have a kind of pattern. This is an indication of further seasonality or periodicity we have not captured. But we can see that the error, trend and seasonality amplitudes does not increase or decrease with time so our standard assumption of additive time series is validated.
We will try to capture the additional seasonality be using the msts function. According to this Wikipedia article there is actually a 22 year magnetic cycle related to the 11 year sunspot cycle. The sunspots are, after all results of the sun’s magnetic activity. This 22 year cycle will be added to the 11 year cycle.
sunspot_11yr_22yr <- msts(data = sun_clean$Monthly.Mean.Total.Sunspot.Number, seasonal.periods = c(12*11, 12*22))
sunspot_11yr_22yr %>% mstl() %>% autoplot()
Now the trend is free of small patterns but it still show cyclical behavior like the sunspot_ts_11yr plot. This is expected as the original plot also shows the count values oscillating around a mean of 0 instead of steadily increasing or decreasing. The remainder does not show regularities like the previous one but there is still a kind of repeated peaks there. There are two likely explanations for the remaining seasonality observed in the remainder component:
The 11-year and 22-year seasonality is not exact but only on average, that is the exact number may vary from 9 to 14 years from decade to decade (link). Additionally the two 11-year cycles within a 22-year cycle also does not repeat exactly.
There are possibilities of longer sunspot cycles (link). Sunspot researchers have not reached any consensus regarding the longer sunspot cycle though.
Another observation is that after the decomposition into 2 seasonalities the 11-year seasonality seems to have 2 different time ranges with one time range have a shorter amplitude than the other.
For the second possibility that there are longer sunspot cycles, let us use the little ice age period, as the frequency. That ice age lasted from 1645 - 1715 that is 70 years.
sunspot_11yr_22yr_70yr <- msts(data = sun_clean$Monthly.Mean.Total.Sunspot.Number, seasonal.periods = c(12*11, 12*22, 12*70))
sunspot_11yr_22yr_70yr %>% mstl() %>% autoplot()
It seems there is a seasonality component that repeats almost perfectly for 70 years. We will use this time series to compare the performance of the baseline forecasting model which is the 11yr_22yr.
4 Cross-Validation
We will separate the last 11 years data for testing and use the rest for training.
sunspot_11yr_22yr_train <- sunspot_11yr_22yr %>% head(-(12*11))
sunspot_11yr_22yr_test <- sunspot_11yr_22yr %>% tail(12*11)5 Forecasting Model
We will use stlm with the arima method.
model_11yr_22yr <- sunspot_11yr_22yr_train %>% stlm(
method = "arima")
mdl_11yr_22yr_f <- forecast(model_11yr_22yr, h=12*11)
accuracy(mdl_11yr_22yr_f$mean, sunspot_11yr_22yr_test)#> ME RMSE MAE MPE MAPE ACF1 Theil's U
#> Test set -13.99055 26.63534 19.30469 -63.50217 158.0615 0.5929517 1.201166autoplot(sunspot_11yr_22yr %>% tail(12*11*5))+
autolayer(mdl_11yr_22yr_f$mean, series = "prediction")+
autolayer(sunspot_11yr_22yr_test, series = "actual test data")
The MAPE is very high, 158% but the autoplot shows that the predicted time series still follow the shape of the actual data although there are some points in which the actual data is low while the predicted data is very high. Furthermore there are some point which are predicted to have negative counts. So although the model cannot predict the number of sunspots with high accuracy, but it can be used for a qualitative prediction of what years during the 11 year cycle will have the maximum sunspot count and minimum count.
6 Attempts to improve the prediction
6.1 Adding new seasonality
If we add a seasonality of 70 years, will the accuracy improve or get worse ?
sunspot_11yr_22yr_70yr_train <- sunspot_11yr_22yr_70yr %>% head(-(11*12))
sunspot_11yr_22yr_70yr_test <- sunspot_11yr_22yr_70yr %>% tail(11*12)
model_11yr_22yr_70yr <- sunspot_11yr_22yr_70yr_train %>% stlm(
method = "arima")
mdl_11yr_22yr_70yr_f <- forecast(model_11yr_22yr_70yr, h=12*11)
accuracy(mdl_11yr_22yr_70yr_f$mean, sunspot_11yr_22yr_70yr_test)#> ME RMSE MAE MPE MAPE ACF1 Theil's U
#> Test set -31.43851 45.44477 36.29108 -169.0625 274.7478 0.701045 2.468598autoplot(sunspot_11yr_22yr_70yr %>% tail(12*11*5))+
autolayer(mdl_11yr_22yr_70yr_f$mean, series = "prediction")+
autolayer(sunspot_11yr_22yr_70yr_test, series = "actual test data")
The MAPE is even worse at 274%. Although once again the general shape is correct, the prediction generally overcounts the sunspot compared to test data. Near the very end, once again the model predicted negative counts. So the model with 2 seasonalities is the better model and adding seasonality does not improve the model.
6.2 Subsetting from the whole data
Next, we will try to subset the train-test data from the whole sunspot data. We notice that the 11 years seasonality and 22 years seasonality have 2 distinct amplitudes and patterns. For the 11 years seasonality the amplitude changes at the sixth peak while for the 22 year seasonality the pattern changes at the eleventh peak from the head or fourteenth peak counted from the tail. So we will use the last 14 peaks and forecast the last peak using the 13 other peaks.
newdata <- sunspot_11yr_22yr %>% tail(14*11*12)
autoplot(mstl(newdata))
Now we can see that the 11 years seasonality has constant amplitudes and the 22 years seasonality has constant patterns.
newdata_train <- newdata %>% head(-(11*12))
newdata_test <- newdata %>% tail(11*12)
model2_11yr_22yr <- newdata_train %>% stlm(
method = "arima")
mdl2_11yr_22yr_f <- forecast(model2_11yr_22yr, h=11*12)
accuracy(mdl2_11yr_22yr_f$mean, newdata_test)#> ME RMSE MAE MPE MAPE ACF1 Theil's U
#> Test set -13.56006 26.22463 19.24945 -103.9421 190.2234 0.5772898 1.695063autoplot(newdata %>% tail(12*11*5))+
autolayer(mdl2_11yr_22yr_f$mean, series = "prediction")+
autolayer(newdata_test, series = "actual test data")
The MAPE is higher than the baseline model at 190% compared to 158%. Once again it follows the general shape of the counts.
6.3 Trying to predict a shorter time period
Using the baseline model, let us try to predict the monthly mean sunspot count for 1 year only
sunspot_11yr_22yr_train2 <- sunspot_11yr_22yr %>% head(-(12))
sunspot_11yr_22yr_test2 <- sunspot_11yr_22yr %>% tail(12)
model3_11yr_22yr <- sunspot_11yr_22yr_train2 %>% stlm(
method = "arima")
mdl3_11yr_22yr_f <- forecast(model3_11yr_22yr, h=12*11)
accuracy(mdl_11yr_22yr_f$mean, sunspot_11yr_22yr_test2)#> ME RMSE MAE MPE MAPE ACF1 Theil's U
#> Test set -0.7941978 10.17434 7.104288 -329.0216 684.2383 0.1767837 1.05579The MAPE becomes very high, 684% but the RMSE decreases from 26 counts of the baseline to 10 counts. This means that If I use this model to predict, then the actual value can be plus-minus 10 counts compared to 26 counts using the baseline model. This is an interesting observation but because we have decided to use MAPE as our metric, this model has to be rejected in favor of the baseline model.
7 Assumptions check
Now let us check the assumptions of the time series best model, there are no autocorrelations between the residuals and that the residuals are distributed normally.
Autocorrelation assumption:
- H_0: The residuals of the model are not autocorrelated
- H_1: The residuals of the model are autocorrelated
Normal Distribution assumption (normality):
- H_0: The residuals of the model are normally distributed
- H_1: The residuals of the model are not normally distributed.
Box.test(model_11yr_22yr$residuals, type = "Ljung-Box")#>
#> Box-Ljung test
#>
#> data: model_11yr_22yr$residuals
#> X-squared = 0.0076018, df = 1, p-value = 0.9305shapiro.test(model_11yr_22yr$residuals)#>
#> Shapiro-Wilk normality test
#>
#> data: model_11yr_22yr$residuals
#> W = 0.97218, p-value < 0.00000000000000022hist(model_11yr_22yr$residuals, breaks=50)
According to the Ljung-Box test with p-value > 0.05 the no autocorrelation of residuals assumption is valid. But, according to the shapiro-wilk test the distribution of the residuals is not normal. If we visualise the residual distribution using histogram, we find that the distribution looks normal. This contrast between low p-value of shapiro-wilk and the bell-curve form of the histogram cannot currently be explained. The Shapiro-Wilk test can take up to 5000 samples, the length of residuals of our model is slightly more than 3000.
length(model_11yr_22yr$residuals)#> [1] 3133So, the shapiro-wilk test should be able to handle the model. In this case we will take the residual histogram over the shapiro-wilk statistical test and conclude that our model satisfied the normality assumption too.
8 Conclusion
We attempted to use time series analysis and forecasting to forecast the number of average monthly number of sunspots for a sunspot cycle, that is, 11 years. Our best model is a multi-seasonal model with arima method and seasonality of 11 and 22 years. It is able to predict the shape of the rise and fall of the sunspot counts during the 11 year period but with an error of more than 100%. Our model then is good enough to make a qualitative prediction but not a quantitative one. We also tried to improve the model by adding one more seasonality of 70 years, subsetting the data before cross validation in order to remove the data with slightly different seasonality pattern and amplitude, and shortening the forecasting period from 11 years to 1 year but the MAPE did not decrease.