# Getting Started ####
library("tidyverse")
library(GGally)
library(ggpubr)
library(ggplot2)
library(olsrr)
library(lm.beta)
library("mctest")
library(splines)
library(forecast)
library("brms")
options(scipen=100, digits=12)
options(max.print = 500000000)
# Simulating Data ####
set.seed(25)
Budget<-runif(1000,10000,25000)
Spend_1<-runif(1000,.10*Budget,.30*Budget)
Spend_2<-runif(1000,.10*Budget,.20*Budget)
Spend_3<-runif(1000,.10*Budget,.40*Budget)
Spend_4<-runif(1000,.01*Budget,.10*Budget)
Actual_Spend<-Spend_1+Spend_2+Spend_3+Spend_4
Week<-1:1000
Sales<-rgamma(1000, shape = 2, rate = 1)*(.004*Actual_Spend)
cor.test(Actual_Spend,Sales)
DF<-data.frame(Budget,Spend_1,Spend_2,Spend_3,Spend_4,Actual_Spend,Sales,Week)
ggpairs(DF)+ggtitle("")+
theme(legend.key = element_blank(),
line = element_line(colour = "black", size = 1),
axis.line = element_line(colour = "black"),
panel.border= element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
title = element_text(family="Times New Roman", size= 24),
axis.title.x = element_text(family="Times New Roman",colour = "Black", margin = margin(t = 20, r = 0, b = 0, l = 0), size = 12),
axis.title.y = element_text(family="Times New Roman",colour = "Black", margin = margin(t = 0, r = 20, b = 0, l = 0), size = 12),
axis.text= element_text(family="Times New Roman", size= 12, color = "Black"),
text= element_text(family="Times New Roman", size= 12),
plot.margin = margin(t = 2, r = 2, b = 2, l = 2, unit = "cm"))
#__________________________________________________________####
# Bayesian Time Series Model ####
#__________________________________________________________####
# Get Priors ####
get_prior(Sales ~ 0 + Intercept+Spend_1+Spend_2+Spend_3+Spend_4, data = DF)
# Specify Priors ####
Priors<-c(set_prior("normal(0, 2.5)", class = "b", coef = "Spend_1"),
set_prior("normal(0, 2.5)", class = "b", coef = "Spend_2"),
set_prior("normal(0, 2.5)", class = "b", coef = "Spend_3"),
set_prior("normal(0, 2.5)", class = "b", coef = "Spend_4"),
set_prior("normal(0, 2.5)", class = "b", coef = "Week"),
set_prior("student_t(3, 0, 56.4)", class = "b", coef = "Intercept"))
M1<-brm(formula =Sales ~ 0 + Intercept+Spend_1+Spend_2+Spend_3+Spend_4+Week, data = DF,
prior = Priors, iter = 5000, warmup = 2000, chains = 4, cores = 12,
control = list(adapt_delta = 0.95, max_treedepth = 15))
print(M1, digits = 10)
print(summary(M1), digits = 5)
#"Odds Ratios"
print(exp(fixef(M1)), digits = 5)
#"model R-squared"
bayes_R2(M1)
rsq_M4<-bayes_R2(M1)
#"Median R-squared for model"
print(median(rsq_M1), digit = 5) # Median rsq for model
#"summary of residuals"
summary(residuals(M1))
#"Correlation matrix to test multicolinearity"
cov2cor(vcov(M1))
#"LOO"
loo(M1)
#"waic"
waic(M1)
#Plot Check
pairs(M1)
#MCMC trace
mcmc_trace(M1)
#Posterior Predictive Check
pp_check(M1)
# Posterior Predict ####
New_Data<-with(DF, data.frame(Week=1001:1052,
Budget = rep(0,52),
Spend_1 = runif(52,30000,35000),
Spend_2 = runif(52,10000,15000),
Spend_3 = rep(0,52),
Spend_4 = runif(52,1000,3000),
Sales = NA))
BP_DF<-data.frame(New_Data, t(posterior_predict(M1, newdata=New_Data, seed = 265810, ndraws = 1000)))
names(BP_DF)
names(BP_DF[8:1007])
BP_DF$Sales<-rowMeans(BP_DF[,c(8:1007)])
row_sd <- apply(BP_DF[,c(8:1007)], 1, sd)
BP_DF$Sales_LB<-BP_DF$Sales-2.60*row_sd
BP_DF$Sales_UB<-BP_DF$Sales+2.60*row_sd
DF1<-(BP_DF[, c("Week","Budget", "Spend_1", "Spend_2","Spend_3","Spend_4","Sales","Sales_LB","Sales_UB")])
DF2<-(DF[, c("Week","Budget", "Spend_1", "Spend_2","Spend_3","Spend_4","Sales")])
DF2$Sales_LB<-NA
DF2$Sales_UB<-NA
Forecast_Spend <- rbind(DF1,DF2)
ggplot(data = Forecast_Spend, aes(x = Week, y = Sales)) +
geom_line(color = 'black', size = 0.7) +
geom_ribbon(aes(ymin = Sales_LB, ymax = Sales_UB), fill = 'red', alpha = 0.2) +
geom_line(color = 'red', size = 0.7, linetype = 'dashed') +
labs(title = "Forecasted") +
theme_minimal()
ianasilver.com 