Forecasting Bitcoin Search Interest

library(tidyverse)
library(lubridate)
library(forecast)
library(tseries)
library(knitr)
library(htmltools)

knitr::opts_chunk$set(
  message = FALSE,
  warning = FALSE,
  fig.align = "center",
  fig.width = 9.5,
  fig.height = 5.2
)

plot_theme <- theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold", size = 16, color = "#1d3557"),
    plot.subtitle = element_text(size = 11, color = "#6c7a89"),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_line(color = "#e6edf5"),
    axis.title = element_text(color = "#334e68"),
    axis.text = element_text(color = "#4b6478")
  )

1. Introduction

This project analyses monthly Google search interest for Bitcoin as a time series. The aim is to investigate whether public interest in Bitcoin contains enough structure to support short-term forecasting.

Bitcoin search interest is not likely to behave like a stable economic or demographic series. Public attention can change quickly when Bitcoin appears in the news, when prices rise or fall sharply, or when cryptocurrency becomes part of wider market discussion. Because of this, the project is not only about producing a forecast, but also about testing how reliable that forecast is.

The analysis uses decomposition, stationarity testing, benchmark models, ETS models and SARIMA models. Models are compared using a train-test split, so the final judgement is based on performance on unseen data rather than visual fit.

Project focus: The main question is whether Bitcoin search interest behaves like a forecastable time series, or whether it is too event-driven for classical models to predict reliably.

2. Data

The dataset contains monthly Bitcoin search interest values. The values are relative Google Trends scores, so they do not represent exact search volumes. A higher value means Bitcoin had higher search interest during that month.

Code: load the Bitcoin search interest dataset

bitcoin_raw <- read.csv("BitcoinTimeSeries.csv")

The original columns are renamed to make the analysis easier to read. The date column is also converted into a proper date format so that R can treat it as a monthly time series.

Code: clean the date column and prepare the search interest series

bitcoin <- bitcoin_raw %>%
  rename(
    date = Time,
    value = Bitcoin
  )

bitcoin$date <- as.Date(bitcoin$date, format = "%d/%m/%Y")

bitcoin <- bitcoin %>%
  filter(!is.na(date), !is.na(value)) %>%
  arrange(date)

start_date <- min(bitcoin$date)
end_date <- max(bitcoin$date)
n_obs <- nrow(bitcoin)
min_value <- min(bitcoin$value)
max_value <- max(bitcoin$value)

After cleaning, the dataset runs from 2015-01-01 to 2026-01-01, with 133 monthly observations. Search interest ranges from 2 to 100.

The data is monthly, so it is converted into a time series object with frequency 12. This tells R that the data has 12 observations per year.

Code: convert the cleaned data into a monthly time series object

btc_ts <- ts(
  bitcoin$value,
  frequency = 12,
  start = c(year(min(bitcoin$date)), month(min(bitcoin$date)))
)

ts_frequency <- frequency(btc_ts)
ts_observations <- length(btc_ts)
ts_start <- paste(start(btc_ts), collapse = "-")
ts_end <- paste(end(btc_ts), collapse = "-")

Frequency

12

Observations

133

Start

2015-1

End

2026-1

3. Exploratory analysis

The first step is to inspect the overall shape of the series.

Code: plot monthly Bitcoin search interest over time

ggplot(bitcoin, aes(x = date, y = value)) +
  geom_line(linewidth = 0.8, colour = "#1d3557") +
  labs(
    title = "Bitcoin search interest over time",
    subtitle = "Monthly relative Google Trends score",
    x = "Date",
    y = "Search interest"
  ) +
  plot_theme +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

The plot shows that search interest is not stable over time. There are several sharp spikes, which is expected for Bitcoin because public attention often increases suddenly around price movements, news events or wider market discussion.

Code: calculate summary statistics for the search interest series

summary_stats <- bitcoin %>%
  summarise(
    observations = n(),
    mean_value = round(mean(value), 2),
    median_value = round(median(value), 2),
    sd_value = round(sd(value), 2),
    minimum_value = min(value),
    maximum_value = max(value)
  )

obs_value <- summary_stats$observations
mean_value <- summary_stats$mean_value
median_value <- summary_stats$median_value
sd_value <- summary_stats$sd_value
minimum_value <- summary_stats$minimum_value
maximum_value <- summary_stats$maximum_value

Observations

133

Mean

20.48

Median

18

Standard deviation

15.5

Minimum

2

Maximum

100

The gap between the median and the maximum suggests that the series is not evenly spread around a typical level. Instead, there are occasional periods where search interest becomes much higher than usual.

Initial observation: Bitcoin search interest appears volatile and spike-driven. A model that assumes a smooth stable pattern is unlikely to be enough by itself.

4. Decomposition

STL decomposition separates the series into trend, seasonal and remainder components. This helps show whether the movement in the series is mainly long-term trend, repeated seasonality or irregular shocks.

Code: decompose the series into trend, seasonality and remainder

decomp <- stl(btc_ts, s.window = "periodic")
plot(decomp)

The trend component shows longer-term movement in Bitcoin interest. The seasonal component shows repeated monthly effects. The remainder shows irregular movements that are not explained by the trend or seasonality.

For this dataset, the remainder is important because Bitcoin search interest is strongly affected by events. This means that even if there is some trend or seasonality, sudden shocks can still dominate the series.

Method note: Decomposition is not a forecasting model by itself. It is used to understand the structure of the series before fitting models.

5. Stationarity

ARIMA-type models usually work better when the series is stationary. A stationary series has statistical properties that are reasonably stable over time. The original plot suggests that Bitcoin search interest may not be stationary, so the Augmented Dickey-Fuller test is used.

Code: run ADF stationarity tests before and after differencing

adf_original <- adf.test(btc_ts)

diff_btc <- diff(btc_ts)
adf_differenced <- adf.test(diff_btc)

adf_table <- data.frame(
  Series = c("Original series", "Differenced series"),
  `ADF statistic` = c(
    round(as.numeric(adf_original$statistic), 4),
    round(as.numeric(adf_differenced$statistic), 4)
  ),
  `p-value` = c(
    round(adf_original$p.value, 4),
    round(adf_differenced$p.value, 4)
  )
)

kable(adf_table, caption = "ADF stationarity tests")

ADF stationarity tests

Series	ADF.statistic	p.value
Original series	-3.3118	0.0723
Differenced series	-6.0517	0.0100

A lower p-value gives stronger evidence against a unit root. The differenced series is included because differencing is a common way to reduce changes in level and make a series more suitable for ARIMA modelling.

Code: plot the differenced series

autoplot(diff_btc) +
  labs(
    title = "Differenced Bitcoin search interest",
    x = "Time",
    y = "Change in search interest"
  ) +
  plot_theme

Modelling point: Differencing changes the problem slightly. Instead of only modelling the level of search interest, the model also focuses on how search interest changes from one month to the next.

6. Autocorrelation

The ACF and PACF plots help assess whether past values are related to current values. This matters because time series models rely on the idea that the past contains information about the future.

Code: plot ACF and PACF for the differenced series

acf(diff_btc, main = "ACF of differenced series")

pacf(diff_btc, main = "PACF of differenced series")

The ACF shows correlation with previous lags, while the PACF shows the correlation at each lag after accounting for shorter lags. These plots are not used to manually choose the final model, but they are useful as a check before using automatic ARIMA selection.

7. Train-test split

To avoid judging models only on how well they fit the past, the final 12 observations are held back as a test set.

Code: split the series into training and test sets

n <- length(btc_ts)

train <- window(btc_ts, end = time(btc_ts)[n - 12])
test <- window(btc_ts, start = time(btc_ts)[n - 11])

train_obs <- length(train)
test_obs <- length(test)

Training set

121

Test set

12

Why this matters: A train-test split gives a more honest comparison. A model can fit historical data well but still perform badly on future observations.

8. Benchmark models

Before fitting ETS and SARIMA models, simple benchmark models are used. These are important because a more complicated model is not automatically better. It should improve on simple alternatives.

The mean forecast is a simple long-term average. The naive forecast assumes the latest observation is the best guide to the near future. The seasonal naive forecast tests whether the same month in the previous year is useful for prediction.

Code: fit mean, naive and seasonal naive benchmark forecasts

mean_fit <- meanf(train, h = length(test))
naive_fit <- naive(train, h = length(test))
snaive_fit <- snaive(train, h = length(test))

benchmark_table <- data.frame(
  Model = c("Mean forecast", "Naive forecast", "Seasonal naive forecast"),
  Idea = c(
    "Future values are forecast using the training-set average.",
    "Future values are forecast using the most recent observation.",
    "Future values repeat the value from the same month in the previous year."
  )
)

kable(benchmark_table, caption = "Benchmark models")

Benchmark models

Model	Idea
Mean forecast	Future values are forecast using the training-set average.
Naive forecast	Future values are forecast using the most recent observation.
Seasonal naive forecast	Future values repeat the value from the same month in the previous year.

Code: plot benchmark forecasts over the test horizon

autoplot(mean_fit) +
  labs(title = "Mean forecast", x = "Time", y = "Search interest") +
  plot_theme

autoplot(naive_fit) +
  labs(title = "Naive forecast", x = "Time", y = "Search interest") +
  plot_theme

autoplot(snaive_fit) +
  labs(title = "Seasonal naive forecast", x = "Time", y = "Search interest") +
  plot_theme

Benchmark purpose: If ETS or SARIMA cannot beat these simple models on unseen data, then the extra complexity is not adding much value.

9. ETS model

ETS model: ETS stands for Error, Trend and Seasonality. It is based on exponential smoothing, meaning it gives more weight to recent observations while still using the history of the series.

ETS models are useful when the level, trend or seasonal pattern of a series changes over time. In this project, ETS is worth testing because Bitcoin search interest does not appear stable across the full period.

Code: fit the ETS model and generate forecasts

ets_fit <- ets(train)
ets_forecast <- forecast(ets_fit, h = length(test))

Code: create a clean ETS model summary table

ets_table <- data.frame(
  Model = as.character(ets_fit$method),
  AIC = round(ets_fit$aic, 2),
  AICc = round(ets_fit$aicc, 2),
  BIC = round(ets_fit$bic, 2)
)

kable(ets_table, caption = "ETS model summary")

ETS model summary

Model	AIC	AICc	BIC
ETS(M,A,M)	939.71	945.65	987.24

Code: plot ETS forecast against the test period

autoplot(ets_forecast) +
  labs(
    title = "ETS forecast",
    subtitle = "Forecast over the 12-month test period",
    x = "Time",
    y = "Search interest"
  ) +
  plot_theme

The forecast line gives the model’s central estimate. The shaded intervals show uncertainty. Wider intervals mean the model is less confident about future values.

10. SARIMA model

SARIMA model: SARIMA is a seasonal version of ARIMA, which stands for AutoRegressive Integrated Moving Average.

The autoregressive part uses past values of the series. The integrated part means the model can difference the series to handle non-stationarity. The moving-average part uses past forecast errors. The seasonal part allows the model to capture patterns that repeat every 12 months.

Code: fit the SARIMA model and generate forecasts

arima_fit <- auto.arima(train, seasonal = TRUE)
arima_forecast <- forecast(arima_fit, h = length(test))

Code: create a clean SARIMA model summary table

sarima_table <- data.frame(
  Model = as.character(arima_fit),
  AIC = round(AIC(arima_fit), 2),
  AICc = round(arima_fit$aicc, 2),
  BIC = round(BIC(arima_fit), 2)
)

kable(sarima_table, caption = "SARIMA model summary")

SARIMA model summary

Model	AIC	AICc	BIC
ARIMA(0,1,0)	884.52	884.55	887.3

Code: plot SARIMA forecast against the test period

autoplot(arima_forecast) +
  labs(
    title = "SARIMA forecast",
    subtitle = "Forecast over the 12-month test period",
    x = "Time",
    y = "Search interest"
  ) +
  plot_theme

A SARIMA model written as ARIMA(p,d,q)(P,D,Q)[12] means that p, d and q are the non-seasonal terms, while P, D and Q are the seasonal terms. The [12] appears because the data is monthly.

11. Forecast comparison

The models are compared on the test set using RMSE, MAE and MAPE.

RMSE penalises large errors more heavily. MAE gives the average absolute error. MAPE expresses the error as a percentage.

Code: calculate out-of-sample accuracy for all models

mean_accuracy <- accuracy(mean_fit, test)
naive_accuracy <- accuracy(naive_fit, test)
snaive_accuracy <- accuracy(snaive_fit, test)
ets_accuracy <- accuracy(ets_forecast, test)
arima_accuracy <- accuracy(arima_forecast, test)

accuracy_table <- rbind(
  Mean = mean_accuracy[2, c("RMSE", "MAE", "MAPE")],
  Naive = naive_accuracy[2, c("RMSE", "MAE", "MAPE")],
  Seasonal_Naive = snaive_accuracy[2, c("RMSE", "MAE", "MAPE")],
  ETS = ets_accuracy[2, c("RMSE", "MAE", "MAPE")],
  SARIMA = arima_accuracy[2, c("RMSE", "MAE", "MAPE")]
)

accuracy_table <- round(as.data.frame(accuracy_table), 3)

kable(accuracy_table, caption = "Out-of-sample forecast accuracy")

Out-of-sample forecast accuracy

	RMSE	MAE	MAPE
Mean	8.686	7.214	24.474
Naive	7.159	6.250	26.222
Seasonal_Naive	8.963	7.333	26.027
ETS	8.103	6.521	24.443
SARIMA	7.159	6.250	26.222

Code: identify the best model using RMSE

best_model_name <- rownames(accuracy_table)[which.min(accuracy_table$RMSE)]
best_rmse <- min(accuracy_table$RMSE)

Best model by RMSE

Naive

Lowest RMSE

7.159

Model selection: The best model is selected using out-of-sample performance, not visual fit. This tests whether the model can generalise beyond the data it was trained on.

12. Residual diagnostics

A model should not leave behind obvious patterns in the residuals. If residuals still show autocorrelation, then the model has not captured all of the time series structure.

Code: run Ljung-Box residual autocorrelation tests

arima_lb <- Box.test(residuals(arima_fit), lag = 12, type = "Ljung-Box")
ets_lb <- Box.test(residuals(ets_fit), lag = 12, type = "Ljung-Box")

residual_table <- data.frame(
  Model = c("ETS", "SARIMA"),
  `Ljung-Box statistic` = c(
    round(as.numeric(ets_lb$statistic), 3),
    round(as.numeric(arima_lb$statistic), 3)
  ),
  `p-value` = c(
    round(ets_lb$p.value, 4),
    round(arima_lb$p.value, 4)
  )
)

kable(residual_table, caption = "Residual autocorrelation diagnostics")

Residual autocorrelation diagnostics

Model	Ljung.Box.statistic	p.value
ETS	6.705	0.8765
SARIMA	8.734	0.7255

The Ljung-Box test checks whether residual autocorrelation remains. A high p-value suggests there is no strong evidence of remaining autocorrelation. A low p-value suggests the model may still be missing some structure.

Code: plot SARIMA residual diagnostics

checkresiduals(arima_fit)

Code: plot ETS residual diagnostics

checkresiduals(ets_fit)

13. Final forecast

After comparing the models, the SARIMA model is refitted on the full dataset and used to produce a 12-month forecast. The train-test split was used for model evaluation; once that is complete, the final model can use all available observations.

Code: refit SARIMA using the full dataset and forecast 12 months ahead

final_model <- auto.arima(btc_ts, seasonal = TRUE)
final_forecast <- forecast(final_model, h = 12)

Code: create final model summary table

final_table <- data.frame(
  Model = as.character(final_model),
  AIC = round(AIC(final_model), 2),
  AICc = round(final_model$aicc, 2),
  BIC = round(BIC(final_model), 2)
)

kable(final_table, caption = "Final model fitted on full dataset")

Final model fitted on full dataset

Model	AIC	AICc	BIC
ARIMA(0,1,0)	964.95	964.98	967.83

Code: plot the final 12-month forecast

autoplot(final_forecast) +
  labs(
    title = "Final 12-month forecast for Bitcoin search interest",
    subtitle = "SARIMA model refitted using the full dataset",
    x = "Time",
    y = "Search interest"
  ) +
  plot_theme

Forecast interpretation: This forecast should not be treated as a precise prediction. Bitcoin search interest is strongly affected by news, price movements and market sentiment. The forecast is best understood as a model-based estimate assuming future behaviour is similar to past patterns.

14. Limitations and possible extensions

There are several limitations to this project.

First, Google Trends gives relative search interest rather than the actual number of searches. Second, this is a univariate analysis, so the models only use past search interest. They do not include Bitcoin price, trading volume, market volatility, news sentiment or macroeconomic variables.

This matters because Bitcoin attention is likely driven by external factors. A large price movement or major news story could change search interest very quickly, and a model based only on past search interest may not react until after the movement has already happened.

Possible extensions include:

• adding Bitcoin price as an external regressor;
• testing an ARIMAX model;
• comparing the results with Prophet;
• using rolling-window cross-validation;
• including news or social media sentiment;
• modelling search interest and price together.

15. Conclusion

This project shows that Bitcoin search interest is volatile, uneven and strongly affected by spikes in public attention. The decomposition and stationarity analysis suggest that the series is not a simple stable process.

The model comparison shows why benchmarks are important. A more complicated model is only useful if it improves on simple methods when tested on unseen data. ETS and SARIMA provide more structured approaches, but their usefulness depends on whether they actually reduce forecast error on the test set.

Overall, classical time series models can capture some of the structure in Bitcoin search interest, but there are clear limits. The series is partly predictable from its own history, but it is also heavily influenced by external events. A stronger future version of this project would include explanatory variables such as Bitcoin price, trading volume and news sentiment.