## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggplot2)
library(tseries)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
##
## Attaching package: 'seasonal'
##
## The following object is masked from 'package:tibble':
##
## view
library(forecast)
library(rmarkdown)
library(knitr)
# Load and inspect the data
data <- read.csv('/Users/kimshellenberger/Desktop/medical_time_series .csv')
head(data)
## Day Revenue
## 1 1 0.0000000
## 2 2 -0.2923555
## 3 3 -0.3277718
## 4 4 -0.3399871
## 5 5 -0.1248875
## 6 6 -0.4915896
## 'data.frame': 731 obs. of 2 variables:
## $ Day : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Revenue: num 0 -0.292 -0.328 -0.34 -0.125 ...
###III: Data Preparation
#C1.Provide a line graph visualizing the realization of the time series.
# 1. Line graph visualization
ggplot(data, aes(x = Day, y = Revenue)) +
geom_line() +
labs(title = "Time Series Visualization",
x = "Day",
y = "Revenue")
####2.Describe the time step formatting of the realization, including any gaps in measurement and the length of the sequence.
#The time step formatting is in days which equals 2 years of data
# Check for any gaps in measurement
print(paste("Any gaps in measurement:", any(diff(data$Day) != time_step)))
## Error in eval(expr, envir, enclos): object 'time_step' not found
# Length of the sequence
length_sequence <- nrow(data)
# 3.Evaluate the stationarity of the time series.
adf.test(data$Revenue)
##
## Augmented Dickey-Fuller Test
##
## data: data$Revenue
## Dickey-Fuller = -2.5906, Lag order = 9, p-value = 0.3283
## alternative hypothesis: stationary
# Take first-order difference to remove trend
data$diff_revenue <- c(NA, diff(data$Revenue))
# Evaluate stationarity
adf_result <- adf.test(data$diff_revenue[-1], alternative = "stationary")
## Warning in adf.test(data$diff_revenue[-1], alternative = "stationary"): p-value
## smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: data$diff_revenue[-1]
## Dickey-Fuller = -7.2794, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary
###The p-value is below .05 so that indicates stationary
# 4.Explain the steps you used to prepare the data for analysis, including the training and test set split.
length_sequence <- nrow(data)
train_data <- data[1:700,]
test_data <- data[-(1:700),]
summary(test_data)
## Day Revenue diff_revenue
## Min. :701.0 Min. :15.25 Min. :-1.17889
## 1st Qu.:708.5 1st Qu.:16.03 1st Qu.:-0.27805
## Median :716.0 Median :18.51 Median : 0.09275
## Mean :716.0 Mean :17.92 Mean :-0.02011
## 3rd Qu.:723.5 3rd Qu.:19.03 3rd Qu.: 0.30914
## Max. :731.0 Max. :20.31 Max. : 1.14450
####I checked for null and n/a in the original data prior to splitting the data into training and test.
#See the 2 attached files.
#The training and test data is split 80/20
# Write train data to CSV
write.csv(train_data, file = "/Users/kimshellenberger/Desktop/train_data.csv", row.names = FALSE)
# Write test data to CSV
write.csv(test_data, file = "/Users/kimshellenberger/Desktop/test_data.csv", row.names = FALSE)
# Part IV: Model Identification and Analysis
# D.1.Report the annotated findings with visualizations of your data analysis, including the following elements:
# •the decomposed time series •the presence or lack of a seasonal component •trends
decomposed_ts <- decompose(ts(data$Revenue, frequency = 2), "multiplicative")
# Plot the decomposed components
plot(decomposed_ts)
# 2.Identify an autoregressive integrated moving average (ARIMA) model that accounts for the observed trend and seasonality of the time series data.
arima_model <- auto.arima(train_data$Revenue)
summary(arima_model)
## Series: train_data$Revenue
## ARIMA(1,1,0)
##
## Coefficients:
## ar1
## 0.4109
## s.e. 0.0345
##
## sigma^2 = 0.1938: log likelihood = -417.86
## AIC=839.72 AICc=839.74 BIC=848.82
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0144574 0.4395562 0.3543905 -0.2323458 7.729072 0.9246115
## ACF1
## Training set -0.001796746
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0144574 0.4395562 0.3543905 -0.2323458 7.729072 0.9246115
## ACF1
## Training set -0.001796746
plot(arima_model$residuals)
# Part V: Data Summary and Implications
# E.Summarize your findings and assumptions by doing the following:
# 1. Discuss the results of your data analysis, including the following points:
# •the selection of an ARIMA model
# the ARIMA model was selected using auto.arima and it mad sense due to no seasonality.
# •the model evaluation procedure and error metric
# I used RMSE as the evaluation metric. The closer the outcome to 1 the better the forcast fit the test data.
rmse <- sqrt(mean((test_data$Revenue - forecast_values)^2))
rmse
# 2. Provide an annotated visualization of the forecast of the final model compared to the test set.
plot(forecast_result, main = "Forecast with 95% Confidence Interval")
lines(train_data$Day, train_data$Revenue, col = "blue") # Train data
lines(test_data$Day, test_data$Revenue, col = "red") # Test data
lines(forecast_result$mean, col = "black") # Forecast
lines(forecast_result$lower[,2], col = "green", lty = 2) # 95% confidence interval lower bound
lines(forecast_result$upper[,2], col = "green", lty = 2) # 95% confidence interval upper bound
legend("topleft", legend = c("Train Data", "Test Data", "Forecast", "95% Confidence Interval"),
col = c("blue", "red", "black", "green"), lty = c(1, 1, 1, 2))
