#2.16

# Load the forecast package
library(forecast)
## Warning: package 'forecast' was built under R version 4.2.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
# Use the built-in sunspots dataset
data(sunspots)

# Extract the first 100 observations of the sunspot series
sunspot_series <- sunspots[1:100]

# Subtract the sample mean from the data to center it around zero
sunspot_centered <- sunspot_series - mean(sunspot_series)

# Plot the sample ACF and PACF up to lag 40
acf(sunspot_centered, lag.max = 40, main = "Sample ACF")

pacf(sunspot_centered, lag.max = 40, main = "Sample PACF")

# Fit an AR(2) model to the centered data
ar2_model <- Arima(sunspot_centered, order = c(2, 0, 0), include.mean = FALSE, method = "ML")

# View the model parameters
print(ar2_model)
## Series: sunspot_centered 
## ARIMA(2,0,0) with zero mean 
## 
## Coefficients:
##          ar1     ar2
##       0.5779  0.3141
## s.e.  0.0938  0.0941
## 
## sigma^2 = 240.2:  log likelihood = -415.68
## AIC=837.37   AICc=837.62   BIC=845.18
# Store the estimated parameters
phi1 <- ar2_model$coef[1]
phi2 <- ar2_model$coef[2]
sigma2 <- ar2_model$sigma2

# Compare the model and sample ACF and PACF
Acf(ar2_model$residuals, main = "Model ACF")

Pacf(ar2_model$residuals, main = "Model PACF")

#2.17

# Use the AR(2) model fitted in the preceding problem
# Forecast the next 10 values
forecast_values <- forecast(ar2_model, h = 10)

# Add the mean to the forecasts
forecast_values_with_mean <- forecast_values$mean + mean(sunspot_series)

# Print the forecast values
print(forecast_values_with_mean)
## Time Series:
## Start = 101 
## End = 110 
## Frequency = 1 
##  [1] 29.84447 30.94805 31.53696 32.22390 32.80584 33.35790 33.85970 34.32307
##  [9] 34.74846 35.13983
# Plot the original data with the forecasts
plot(forecast_values_with_mean, xlab = "Time", ylab = "Sunspot values",
     main = "Sunspot Series with Forecasts", xlim = c(101, 110))
lines(sunspot_series, col = "blue")
legend("topleft", legend = c("Original Data", "Forecasts"), col = c("blue", "red"), lty = 1)