Monthly Boxplot of Bitcoin Prices

The range of prices varies significantly across the months, with some months showing higher variability and higher maximum prices. Month 4 and month 6 show the highest median prices and also a wide range of prices, indicating significant volatility. Month 9 and Month 10 show lower median prices and a relatively narrower price range. Month 2 and Month 3 have moderate median prices with a moderate range of variation. the plot provides a clear visual representation of how Bitcoin prices fluctuate month by month, showing periods of high volatility and identifying months with significant price changes.

Yearly Boxplot of Bitcoin Prices

The yearly boxplot shows the lowest bitcoin prices with minimal fluctuation in the year 2015-2016. 2017-2019 shows a gradual increase in prices with some fluctuations over the months. 2020 is interesting because although the median shows gradual increase but a few points are way above the 75th percentile, indicating he start of high fluctuations. 2021 is the peak with a gradual decline in 2022.

Year-wise Trend Lines of Bitcoin Prices

This plot provides a visual representation of the price trends over time. It helps in identifying overall growth patterns, significant peaks, and drops. Once again, the years 2020-to 2022 show great fluctuations and the peak of bitcoin prices.

Time series of Monthly Bitcoin Prices

The time series clearly shows an upward trend and slight seasonality from the end of the year 2000.

Correlation between Consecutive Months

A high correlation (0.9618) between consecutive months suggests that Bitcoin prices are highly dependent on the prices in the previous month, indicating a strong autocorrelation in the data.

4.3.1 Linear Regression

# Create a linear model of the time series dataset
linear_model <- lm(Price ~ Date, data = BitCoin_df)

# Show the summary of the model and explain the outcome
summary(linear_model)

Call:
lm(formula = Price ~ Date, data = BitCoin_df)

Residuals:
   Min     1Q Median     3Q    Max 
-15114  -7997  -2255   3065  35626 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -2.211e+05  1.939e+04  -11.40   <2e-16 ***
Date         1.308e+01  1.073e+00   12.19   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10430 on 105 degrees of freedom
Multiple R-squared:  0.586, Adjusted R-squared:  0.5821 
F-statistic: 148.6 on 1 and 105 DF,  p-value: < 2.2e-16

Regression Analysis

# Create a plot of the linear model on top of the time series dataset line plot with scatter data points
ggplot(BitCoin_df, aes(x = Date, y = Price)) + 
  geom_point() + 
  geom_line() + 
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(x = "Time", y = "Monthly Average BTC Price") +
  theme_minimal() + 
  ggtitle("Linear Regression Model on Time Series Data")

# Perform residual analysis and create a line & scatter plot of the residuals. Explain the outcome
residuals <- resid(linear_model)
plot(BitCoin_df$Date, residuals,xlab="Time", ylab="Residuals", type = "o", col = "blue", main = "Residuals of Linear Model")
abline(h=0.05, col= "red",lty=2)

# Create a histogram plot of the residuals. Explain the outcome
hist(residuals, xlab="Residuals", breaks = 30, col = "lightblue", main = "Histogram of Residuals")

# Create ACF & PACF plots of residuals. Explain the outcome
acf(residuals)

pacf(residuals)

# Create QQ plot of residuals. Explain the outcome
qqnorm(residuals)
qqline(residuals, col = "red")

# Perform Shapiro-Wilk test on residuals. Explain the outcome
shapiro.test(residuals)

    Shapiro-Wilk normality test

data:  residuals
W = 0.85983, p-value = 1.215e-08

4.3.2 Quadratic Regression

# Create a quadratic model of the time series dataset
BitCoin_df$Date <- as.numeric(BitCoin_df$Date)
quadratic_model <- lm(Price ~ poly(Date, 2), data = BitCoin_df)

# Show the summary of the model and explain the outcome
summary(quadratic_model)

Call:
lm(formula = Price ~ poly(Date, 2), data = BitCoin_df)

Residuals:
   Min     1Q Median     3Q    Max 
-15872  -7420  -1996   2666  36106 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)       14944       1010  14.794   <2e-16 ***
poly(Date, 2)1   127161      10449  12.170   <2e-16 ***
poly(Date, 2)2     8246      10449   0.789    0.432    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10450 on 104 degrees of freedom
Multiple R-squared:  0.5885,    Adjusted R-squared:  0.5806 
F-statistic: 74.36 on 2 and 104 DF,  p-value: < 2.2e-16

Testing The Data Set for Stationarity

# Load necessary libraries
library(lmtest)

# Convert the Bitcoin data frame to a time series object with frequency 12 (monthly data)
btc_ts <- ts(BitCoin$Price, start = c(2015, 1), frequency = 12)

# Check for and handle missing values
if (any(is.na(btc_ts))) {
  btc_ts <- na.approx(btc_ts)  # Linear interpolation to handle missing values
}
# ACF & PACF
par(mfrow=c(1,2))
acf(btc_ts, lag.max = 24)
pacf(btc_ts, lag.max = 24)

• ACF and PACF Explanation: The ACF plot shows a slow decay of autocorrelations, indicating non-stationarity and the need for differencing to achieve stationarity, common in financial time series. The PACF plot shows a significant spike at lag 1, suggesting the series may follow an AR(1) process, indicating that the AR term in the ARIMA model may be sufficient with a low order.

library(tseries)
adf.test(btc_ts)

    Augmented Dickey-Fuller Test

data:  btc_ts
Dickey-Fuller = -2.5743, Lag order = 4, p-value = 0.3385
alternative hypothesis: stationary

• ADF Explanation: The Augmented Dickey-Fuller test results for the Bitcoin time series (btc_ts) indicate non-stationarity (p-value = 0.3385), suggesting the need for differencing to achieve stationarity before fitting an ARIMA model.

qqnorm(y=btc_ts,main = "QQ Plot.")
qqline(y=btc_ts,col=2,lwd=1,lty=2)

- QQPlot Explanation of Bitcoin Time Series: The QQ plot shows that the sample quantiles deviate significantly from the theoretical quantiles, especially in the tails, indicating that the Bitcoin time series data is not normally distributed. The heavy tails suggest the presence of extreme values or outliers. This non-normality needs to be addressed, possibly by transforming the data or using models that do not assume normality.

sw <-  shapiro.test(btc_ts)
sw

    Shapiro-Wilk normality test

data:  btc_ts
W = 0.83358, p-value = 1.258e-09

• Shapiro-Wilk Test Explanation of Bitcoin Time Series:

Since the p-value is very small (much less than the typical significance level of 0.05), we reject the null hypothesis that the data is normally distributed. This suggests that btc_ts does not follow a normal distribution.

Differencing The Data for Stationarity

#Differencing
stationary <- diff(btc_ts)
par(mfrow=c(1,1))
plot(stationary,type='o',ylab="value series",main="Times series plot of the first difference of the Bitcoin price series")

adf.test(stationary)

    Augmented Dickey-Fuller Test

data:  stationary
Dickey-Fuller = -5.1599, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary

#Model Selection
par(mfrow=c(1,2))
acf(stationary, main= "ACF plot of the first difference", lag.max = 30)
pacf(stationary, main= "PACF plot of the first difference", lag.max = 30)

The provided ACF (Autocorrelation Function) and PACF (Partial Autocorrelation Function) plots of the first difference of Bitcoin prices suggest that the time series is now stationary. In the ACF plot, the significant spike at lag 1 and the rapid decline thereafter indicate that the differencing has successfully removed the trend. In the PACF plot, the significant spike at lag 1 followed by smaller spikes suggests that an AR(1) model might be appropriate for the differenced series. Overall, these plots help in identifying the order of ARIMA models for time series forecasting.

#EACF test
library(TSA)
eacf(stationary, ar.max=3, ma.max=3)
AR/MA
  0 1 2 3
0 o o o o
1 x x o o
2 o x o o
3 x o x o

Based on the EACF, we select the ARIMA parameters of models (0,1,1), (0,1,2), (3,1,1)

#ARIMA(0,1,1)
model_011= arima(stationary,order=c(0,1,1))
coeftest(model_011)

z test of coefficients:

     Estimate Std. Error z value  Pr(>|z|)    
ma1 -1.000000   0.028907 -34.594 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#ARIMA(0,1,2)
model_012= arima(stationary,order=c(0,1,2))
coeftest(model_012)

z test of coefficients:

     Estimate Std. Error z value  Pr(>|z|)    
ma1 -0.735344   0.099068 -7.4227 1.148e-13 ***
ma2 -0.264656   0.095704 -2.7654  0.005686 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_311= arima(stationary,order=c(3, 1, 1))
coeftest(model_311)

z test of coefficients:

     Estimate Std. Error  z value Pr(>|z|)    
ar1  0.232756   0.097733   2.3816  0.01724 *  
ar2 -0.190737   0.098518  -1.9361  0.05286 .  
ar3 -0.021621   0.097960  -0.2207  0.82531    
ma1 -0.999999   0.028776 -34.7507  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Extract AIC and BIC values
aic_values <- c(AIC(model_011), AIC(model_012), AIC(model_311))
bic_values <- c(BIC(model_011), BIC(model_012), BIC(model_311))

# Create a data frame
model_table <- data.frame(
  Model = c("model_011", "model_012", "model_311"),
  AIC = aic_values,
  BIC = bic_values
)

# Print the table
print(model_table)
      Model      AIC      BIC
1 model_011 2070.910 2076.218
2 model_012 2066.781 2074.743
3 model_311 2068.434 2081.703

Best Models Based on AIC and BIC:

Model 012 (ARIMA(0,1,2)) is the best choice according to both AIC and BIC, indicating a strong candidate for selection.

Second Model: Model 311 (ARIMA(3,1,1)) has the second lowest AIC and a relatively low BIC.

Therefore, while both AIC and BIC point to Model 012 as the best fit, Model 311 is also a reasonable choice, especially if predictive performance is a higher priority than model simplicity.

model_012 = Arima(stationary, order = c(0, 1, 2))
model_311 = Arima(stationary, order = c(3, 1, 1))

a012 <- accuracy(model_012)
a311 <- accuracy(model_311)



df_models <- data.frame(
rbind(a012, a311)
)
colnames(df_models) <- c("ME", "RMSE", "MAE", "MPE", "MAPE", "MASE", "ACF1")
rownames(df_models) <- c("ARIMA(0,1,2)",  "ARIMA(3,1,1)")

kable(df_models, digits = 2, formats="html", row.names = TRUE) %>%
  kable_styling(full_width = F, font_size = 12, position = "center")

The ARIMA(3,1,1) model has a slightly lower RMSE (4254.90 vs. 4312.32) and MAE (2423.61 vs. 2490.54) compared to the ARIMA(0,1,2) model, indicating better overall accuracy in terms of these metrics. The ACF1 values for both models are close to zero (-0.03 for ARIMA(0,1,2) and -0.01 for ARIMA(3,1,1)), suggesting minimal autocorrelation in the residuals, which is desirable. Despite these differences, both models exhibit high MAPE values, indicating substantial forecasting errors. Overall, ARIMA(3,1,1) performs slightly better in terms of error metrics, but both models have similar performance in terms of residual autocorrelation.

res_012 <- rstandard(model_012)
res_311 <- rstandard(model_311)

#Model 012
par(mfrow = c(1, 1))

# Plot 1: Time series plot of standardized residuals
plot(res_012,
     xlab = "Year Index", ylab = "Residuals",
     main = "Time series plot of standardized residuals-Model 012")
lines(res_012, col = "blue")
abline(h = 0.05, col = "red", lty = 2)

#Model 311
par(mfrow = c(1, 1))

# Plot 1: Time series plot of standardized residuals
plot(res_311,
     xlab = "Year Index", ylab = "Residuals",
     main = "Time series plot of standardized residuals- Model 311")
lines(res_311, col = "blue")
abline(h = 0.05, col = "red", lty = 2)

Both model 012 and 311 show stationarity.

#Histogram 012
hist(res_012,breaks=20, col="skyblue",
     xlab="Residuals",ylab="Frequency",
     main="Histogram of Residuals Model-012")

#Histogram 311
hist(res_311,breaks=20, col="skyblue",
     xlab="Residuals",ylab="Frequency",
     main="Histogram of Residuals Model-311")

Both show normal distribution and very similar spread.

# ACF & PACF of Residuals 012
acf_012 <- acf(res_012)

pacf_012 <- pacf(res_012)

# ACF & PACF of Residuals 311
acf_311 <- acf(res_311)

pacf_311 <- pacf(res_311)

ACF: Both models have residuals with minimal autocorrelation, indicated by most spikes being within the 95% confidence bounds. Given the slightly lower residual variance and error metrics for ARIMA(3,1,1) compared to ARIMA(0,1,2), the ARIMA(3,1,1) model is considered slightly better.

PACF: Both PACF plots indicate minimal partial autocorrelation, with most spikes within the confidence bounds for both models. Given the slightly better performance metrics and lower residual variance observed earlier, the ARIMA(3,1,1) model is considered slightly better overall.

#QQ Plot of Residuals 012
qqnorm(res_012)
qqline(res_012, col="red")

#QQ Plot of Residuals 311
qqnorm(res_311)
qqline(res_311, col="blue")

Both Q-Q plots show that the residuals of the ARIMA(0,1,2) and ARIMA(3,1,1) models are approximately normally distributed with slight deviations at the tails.

# Shapiro-Wilk Test on Residuals 012
print(shapiro.test(res_012))

    Shapiro-Wilk normality test

data:  res_012
W = 0.85275, p-value = 7.244e-09

# Shapiro-Wilk Test on Residuals 311
print(shapiro.test(res_311))

    Shapiro-Wilk normality test

data:  res_311
W = 0.84209, p-value = 2.889e-09

Both models’ residuals fail the Shapiro-Wilk normality test, indicating significant non-normality. The ARIMA(0,1,2) model has a slightly higher W statistic than the ARIMA(3,1,1) model, suggesting that its residuals are closer to normality.

Despite this, when considering previous analyses (lower residual variance and better fit metrics), the ARIMA(3,1,1) model is still considered slightly better overall.

Linear Regression:

Pros: Simple to implement and interpret.

Cons: Limited by linearity assumption, residuals show non-normality and autocorrelation.

Quadratic Regression:

Pros: Can model slight curvature in the trend.

Cons: The quadratic term was not significant, similar performance to linear regression.

ARIMA:

Pros: Captures both trend and seasonality, well-suited for time series data.

Cons: More complex to implement and interpret, requires careful model identification and validation.