Imagine our Sun as a giant cosmic dancer, spinning with dark blemishes that come and go in an 11-year rhythm. These sunspots - temporary dark patches caused by intense magnetic activity - are like the Sun’s heartbeat monitor. Using real solar data, let’s analyze their fascinating patterns.
# Sunspots sunspot_data <-read_delim("https://www.sidc.be/silso/DATA/SN_d_tot_V2.0.csv",delim =";",col_names =c("year", "month", "day", "date_decimal", "ssn", "sd","n_obs", "definitive") )sunspot_data <- sunspot_data %>%# Convert character columns to numeric (keep NAs for filtering)mutate(across(where(is.character), ~parse_number(.) %>%as.numeric())) %>%# Create valid datetime columnmutate(date =make_datetime(year =replace_na(year, 0), # Handle potential NA yearsmonth =replace_na(month, 1), # Default to January if missingday =replace_na(day, 1) # Default to 1st day if missing ),.after = day ) %>%# Filter invalid dates and convert to tsibblefilter(!is.na(date)) %>%as_tsibble(index = date,regular =FALSE# Explicitly declare irregular time series )sunspot_data %>%filter(year(date) >2000) %>%autoplot(ssn) +labs(title ="Daily Sunspot Numbers since 2000",y ="Mean Sunspot Number",x ="Date") +theme_minimal()
The plot shows dramatic spikes and valleys - during solar maximum (2001-2002, 2014-2015) the Sun becomes peppered with spots, while solar minimum (2008-2009, 2019-2020) leaves it nearly spotless.
By averaging to monthly values (since 1800), we see the characteristic sawtooth pattern of solar cycles emerge more clearly. The smoothed curve helps us identify historical extremes like the Dalton Minimum (1790-1830) and Modern Maximum (1950-2000).
# Clean and transform the data firstsolar_analysis <- sunspot_data_monthly %>% tsibble::fill_gaps() %>% tidyr::fill(mean_ssn, .direction ="down")# Create STL decomposition with appropriate parameterssolar_model <- solar_analysis %>%model(STL(mean_ssn ~season(period ="1 year") +trend(),robust =TRUE) )# Visualize componentscomponents(solar_model) %>%autoplot() +labs(title ="Solar Cycle STL Decomposition",subtitle ="11-year Cycle Analysis") +theme_bw()
Solar Cycle Anatomy Lesson
Using STL decomposition, we dissect the signal into three components:
Trend: Gradual intensity changes over centuries
Seasonality: The dominant ~11-year cycle
Remainder: Unexplained “solar tantrums”
# 1. Prepare data with 11-year cycle featuressolar_ts <- sunspot_data_monthly %>%mutate(cycle = ((year(month) -1755) %/%11), # Wolf cycle numberinglog_ssn =log10(mean_ssn +1e-6) ) %>%fill_gaps() %>%fill(log_ssn, .direction ="down") %>%drop_na()# 2. Fit multiple models with Fourier terms for 11-year cyclesolar_models <- solar_ts %>%model(ETS =ETS(log_ssn ~season("A") +trend("N") +error("A")),ARIMA =ARIMA(log_ssn ~pdq() +PDQ() +fourier(period =132, K =5)),STL_ARIMA =decomposition_model(STL(log_ssn ~season(window=99)), ARIMA(season_adjust ~pdq() +fourier(period=132, K=2)) ) )solar_models %>% dplyr::select(ETS) %>%report()
Our forecasting showdown compares three models. The ARIMA model with Fourier terms (MASE=0.578) outperforms others by capturing both cyclical patterns and random fluctuations. All models suggest we’re entering Solar Cycle 25’s ascending phase!
