Saloni-Sarbhai
commented
4 months ago # Load necessary library
library(reprex)
#> Warning: package 'reprex' was built under R version 4.2.3
library(fredr)
#> Warning: package 'fredr' was built under R version 4.2.3
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
library(tidyverse)
#> Warning: package 'tidyverse' was built under R version 4.2.3
#> Warning: package 'tibble' was built under R version 4.2.3
#> Warning: package 'dplyr' was built under R version 4.2.3
#> Warning: package 'forcats' was built under R version 4.2.3
library(lubridate)
library(tseries)
#> Warning: package 'tseries' was built under R version 4.2.3
fredr_set_key("619954596c335c6cd3c4fc1f7a346118")
icnsa <- fredr(series_id = "ICNSA")
# Display the first few rows of the data
head(icnsa_data)
#> Error in head(icnsa_data): object 'icnsa_data' not found
icnsa_data$DATE <- as.Date(icnsa_data$DATE)
#> Error in as.Date(icnsa_data$DATE): object 'icnsa_data' not found
summary(icnsa_data)
#> Error in summary(icnsa_data): object 'icnsa_data' not found
unemployment_data <- fredr(series_id = "UNRATE")
# Display the first few rows of the data
head(unemployment_data)
#> # A tibble: 6 × 5
#> date series_id value realtime_start realtime_end
#> <date> <chr> <dbl> <date> <date>
#> 1 1948-01-01 UNRATE 3.4 2024-02-22 2024-02-22
#> 2 1948-02-01 UNRATE 3.8 2024-02-22 2024-02-22
#> 3 1948-03-01 UNRATE 4 2024-02-22 2024-02-22
#> 4 1948-04-01 UNRATE 3.9 2024-02-22 2024-02-22
#> 5 1948-05-01 UNRATE 3.5 2024-02-22 2024-02-22
#> 6 1948-06-01 UNRATE 3.6 2024-02-22 2024-02-22
unemployment_data$date <- as.Date(unemployment_data$date)
summary(unemployment_data)
#> date series_id value realtime_start
#> Min. :1948-01-01 Length:913 Min. : 2.500 Min. :2024-02-22
#> 1st Qu.:1967-01-01 Class :character 1st Qu.: 4.400 1st Qu.:2024-02-22
#> Median :1986-01-01 Mode :character Median : 5.500 Median :2024-02-22
#> Mean :1985-12-31 Mean : 5.703 Mean :2024-02-22
#> 3rd Qu.:2005-01-01 3rd Qu.: 6.700 3rd Qu.:2024-02-22
#> Max. :2024-01-01 Max. :14.800 Max. :2024-02-22
#> realtime_end
#> Min. :2024-02-22
#> 1st Qu.:2024-02-22
#> Median :2024-02-22
#> Mean :2024-02-22
#> 3rd Qu.:2024-02-22
#> Max. :2024-02-22
library(ggplot2)
# Create the ggplot object
plot <- ggplot(data = icnsa, aes(x = date, y = value)) +
geom_line() +
labs(title = "ICNSA", x = "Year", y = "Number")
# Display the plot
print(plot)
library(ggplot2)
# Create the ggplot object and store it in a variable
plot <- ggplot(data = unemployment_data, aes(x = date, y = value)) +
geom_line() +
labs(title = "Unemployed", x = "Year", y = "Number")
# Display the plot
print(plot)
library(dplyr)
# Mutate to extract Year and Month
icnsa <- icnsa %>%
mutate(Year = lubridate::year(date),
Month = lubridate::month(date))
icnsa_monthly <- icnsa %>%
group_by(Month, Year) %>%
summarise(ICNSA = mean(value, na.rm = TRUE), .groups = "drop")
unemployment_mdata <- unemployment_data %>%
mutate(Year = lubridate::year(date),
Month = lubridate::month(date))
# Merge datasets by Year and Month, retaining all observations
merged_data <- full_join(icnsa_monthly, unemployment_mdata, by = c("Year", "Month"))
# Remove rows with missing values
merged_data <- na.omit(merged_data)
correlation <- cor(merged_data$ICNSA, merged_data$value)
cat("correlation is:",correlation)
#> correlation is: 0.5061069
icsna_ts <- ts(merged_data$ICNSA, frequency = 12, start = c(min(merged_data$Year), 1))
adf_result_icsna <- adf.test(icsna_ts)
#> Warning in adf.test(icsna_ts): p-value smaller than printed p-value
print(adf_result_icsna)
#>
#> Augmented Dickey-Fuller Test
#>
#> data: icsna_ts
#> Dickey-Fuller = -7.8661, Lag order = 8, p-value = 0.01
#> alternative hypothesis: stationary
covariate_ts <- ts(merged_data$value, frequency = 12, start = c(min(merged_data$Year), 1))
# Perform Augmented Dickey-Fuller test for covariate
adf_result_covariate <- adf.test(covariate_ts)
#> Warning in adf.test(covariate_ts): p-value smaller than printed p-value
print(adf_result_covariate)
#>
#> Augmented Dickey-Fuller Test
#>
#> data: covariate_ts
#> Dickey-Fuller = -8.6262, Lag order = 8, p-value = 0.01
#> alternative hypothesis: stationary
final_model <- auto.arima(y = icsna_ts, xreg = covariate_ts)
summary(final_model)
#> Series: icsna_ts
#> Regression with ARIMA(1,1,2) errors
#>
#> Coefficients:
#> ar1 ma1 ma2 xreg
#> 0.6515 -1.3357 0.3446 98138.785
#> s.e. 0.0787 0.0954 0.0930 6912.027
#>
#> sigma^2 = 3.744e+10: log likelihood = -9296.19
#> AIC=18602.39 AICc=18602.48 BIC=18625.03
#>
#> Training set error measures:
#> ME RMSE MAE MPE MAPE MASE
#> Training set -4721.292 192795.5 83964.92 -3.101819 20.8337 0.6004364
#> ACF1
#> Training set -0.0004868408
checkresiduals(final_model)
#>
#> Ljung-Box test
#>
#> data: Residuals from Regression with ARIMA(1,1,2) errors
#> Q* = 25.514, df = 21, p-value = 0.2256
#>
#> Model df: 3. Total lags used: 24
# Generate forecast for the next time point
forecast_values <- forecast(final_model, xreg = covariate_ts, h = 1)
# Extract the mean of the forecasted values for the next time point
single_forecast <- as.numeric(forecast_values$mean[1])
# Print the next predicted value
cat("Next predicted value is", single_forecast, "\n")
#> Next predicted value is 245415.5Created on 2024-02-22 with reprex v2.0.2
Created on 2024-02-09 with reprex v2.0.2