Unemployment Forecast: 3.548072% Method: Bayesian AR6X forecast with external series Initial Unemployment Claims (ISCA) and Industrial Production Index (INDPRO))
Payrolls Forecast: 200.74 Method: AR6x with external series Initial Unemployment Claims (ISCA) and ADP payrolls Survey (NPPTTL)
Motivation
For my final unemployment forecast I chose a Bayesian AR6x model. The reason I chose this model is because it relies heavily on both quantitative and qualitative measures. Prior selection and external series both play a large part in developing a posterior predictive distribution that can accurately predict the unemployment rate. I faced a large series of challenges with the set of Bayesian models I had tested for my final forecast. These included Bayesian Vector Auto regression methods with Normal Wishart,Minnesota and Dirichlet priors. All of which were too computationally intensive to run on kaggle without waiting a very long time for results. The kaggle server even listed that some computations had non-zero exit status. After switching to computing the models in the actual R-Studio software was I even able to run the models efficiently. The selection of priors that I decided to test were based on the literature Volkan Sevin,2009, and Gómez-Rubio,2021. These articles outlined the best priors to use when estimating BVAR for the unemployment rate.
For all of my Bayesian models I used only a portion of the unemployment data past 2002. this was because Bayesian models perform best when there is only a small amount of data and their predictive power becomes less useful when there is a long and extensive history of data. Only testing the data’s rmse against the previous 5 months predictions I was able to achieve an rmse of .26 for my Bayesian AR6x model and all other sets of priors were >.35. Using the Bayesian model really tested my ability to understand what shaping a posterior predictive distribution that explains the data well takes. Experimenting with a variety of different priors was very stressful but also very fulfilling and interesting.
For the payroll forecast I went with a simpler approach. Throughout the semester forecasting the growth in non-farm payroll has been very difficult because of the large variation in the data. With this in mind I wanted to use a simple forecasting method that did not include a variety of complex assumptions and wasn’t too statistically rigorous. The regular ARIMAX model with the initial unemployment claims and ADP payroll Survey seemed to perform well in the past and therefore I thought it would do well for the month of April. My choice of external series was based on tests with a variety of external series used throughout the course.
From Month 1 to April my forecasts moved from simple arima models, the naive forecasts and VARs to Bayesian methods and arima models using external series. As I learned more about forecasting I became more comfortable branching out and trying new models.
library(tidyverse) # metapackage of all tidyverse packages
Warning: package ‘tidyverse’ was built under R version 4.0.5
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
-- Attaching packages ---------------------------------------------------------------------- tidyverse 1.3.1 --
√ ggplot2 3.3.6 √ purrr 0.3.4
√ tibble 3.1.6 √ dplyr 1.0.8
√ tidyr 1.2.0 √ stringr 1.4.0
√ readr 2.1.2 √ forcats 0.5.1
Warning: package ‘tibble’ was built under R version 4.0.5
Warning: package ‘tidyr’ was built under R version 4.0.5
Warning: package ‘readr’ was built under R version 4.0.5
Warning: package ‘dplyr’ was built under R version 4.0.5
Warning: package ‘forcats’ was built under R version 4.0.5
-- Conflicts ------------------------------------------------------------------------- tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
library(fpp2)
Warning: package ‘fpp2’ was built under R version 4.0.5
Registered S3 method overwritten by 'quantmod':
method from
as.zoo.data.frame zoo
-- Attaching packages ----------------------------------------------------------------------------- fpp2 2.4 --
√ forecast 8.16 √ expsmooth 2.3
√ fma 2.4
Warning: package ‘forecast’ was built under R version 4.0.5
Warning: package ‘fma’ was built under R version 4.0.5
Warning: package ‘expsmooth’ was built under R version 4.0.5
library(ggplot2)
library(dynlm)
Warning: package ‘dynlm’ was built under R version 4.0.5
Loading required package: zoo
Warning: package ‘zoo’ was built under R version 4.0.5
Attaching package: ‘zoo’
The following objects are masked from ‘package:base’:
as.Date, as.Date.numeric
library(glmnet)
Warning: package ‘glmnet’ was built under R version 4.0.5
Loading required package: Matrix
Attaching package: ‘Matrix’
The following objects are masked from ‘package:tidyr’:
expand, pack, unpack
Loaded glmnet 4.1-4
library(rstanarm)
Warning: package ‘rstanarm’ was built under R version 4.0.5
Loading required package: Rcpp
Warning: package ‘Rcpp’ was built under R version 4.0.5
This is rstanarm version 2.21.3
- See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!
- Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.
- For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores())
library(fredr)
Warning: package ‘fredr’ was built under R version 4.0.5
library(quantreg)
Warning: package ‘quantreg’ was built under R version 4.0.5
Loading required package: SparseM
Warning: package ‘SparseM’ was built under R version 4.0.4
Attaching package: ‘SparseM’
The following object is masked from ‘package:base’:
backsolve
Warning in .recacheSubclasses(def@className, def, env) :
undefined subclass "packedMatrix" of class "replValueSp"; definition not updated
Warning in .recacheSubclasses(def@className, def, env) :
undefined subclass "packedMatrix" of class "mMatrix"; definition not updated
library(vars)
Warning: package ‘vars’ was built under R version 4.0.5
Loading required package: MASS
Attaching package: ‘MASS’
The following object is masked _by_ ‘.GlobalEnv’:
survey
The following objects are masked from ‘package:fma’:
cement, housing, petrol
The following object is masked from ‘package:dplyr’:
select
Loading required package: strucchange
Warning: package ‘strucchange’ was built under R version 4.0.5
Loading required package: sandwich
Warning: package ‘sandwich’ was built under R version 4.0.5
Attaching package: ‘strucchange’
The following object is masked from ‘package:stringr’:
boundary
Loading required package: urca
Warning: package ‘urca’ was built under R version 4.0.5
Loading required package: lmtest
Warning: package ‘lmtest’ was built under R version 4.0.5
library(prophet)
Warning: package ‘prophet’ was built under R version 4.0.5
Loading required package: rlang
Warning: package ‘rlang’ was built under R version 4.0.5
Attaching package: ‘rlang’
The following objects are masked from ‘package:purrr’:
%@%, as_function, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl, flatten_raw,
invoke, splice
library(xts)
Warning: package ‘xts’ was built under R version 4.0.5
Attaching package: ‘xts’
The following objects are masked from ‘package:dplyr’:
first, last
library(BVAR)
Warning: package ‘BVAR’ was built under R version 4.0.5
Attaching package: ‘BVAR’
The following objects are masked from ‘package:vars’:
fevd, irf
library(Metrics)
Warning: package ‘Metrics’ was built under R version 4.0.5
Attaching package: ‘Metrics’
The following object is masked from ‘package:rlang’:
ll
The following object is masked from ‘package:rstanarm’:
se
The following object is masked from ‘package:forecast’:
accuracy
library(forecast)
library(Boom)
Warning: package ‘Boom’ was built under R version 4.0.5
Attaching package: ‘Boom’
The following object is masked from ‘package:stats’:
rWishart
fredr_set_key("8782f247febb41f291821950cf9118b6")
UNRATE<-fredr(series_id = "UNRATE")
PAYEMS<-fredr(series_id = "PAYEMS",units="chg")
urate<-ts(UNRATE$value,frequency=12,start=c(1948,1),names="Unemployment")
payems<-ts(PAYEMS$value,frequency=12,start=c(1939,1),names="Payrolls")
Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).
The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.
fredr_set_key("8782f247febb41f291821950cf9118b6") #Key I obtained for this class
#The unemployment rate reported in the BLS report has FRED ID "UNRATE"
#Use vintage_dates to get vintage data
UNRATE<-fredr(series_id = "UNRATE",vintage_dates=as.Date("2022-4-25"))
#You can also download this series at https://fred.stlouisfed.org/series/UNRATE
#The jobs series has FRED ID "PAYEMS", for "Payroll Employment Survey"
#Units="chg" ensures we get growth rather than level
#Use vintage_dates to get vintage data (works with changes)
PAYEMS<-fredr(series_id = "PAYEMS",units="chg",vintage_dates=as.Date("2022-4-25"))
#You can also download this series at https://fred.stlouisfed.org/series/PAYEMS
#Format the series as monthly time series objects, starting at the first date
urate<-ts(UNRATE$value,frequency=12,start=c(1948,1),names="Unemployment")
payems<-ts(PAYEMS$value,frequency=12,start=c(1939,1),names="Payrolls")
currentdate<-c(2022,3)
uratelags<-data.frame(
window(urate,start=c(1948,7),end=currentdate),
window(stats::lag(urate,-1),start=c(1948,7),end=currentdate),
window(stats::lag(urate,-2),start=c(1948,7),end=currentdate),
window(stats::lag(urate,-3),start=c(1948,7),end=currentdate),
window(stats::lag(urate,-4),start=c(1948,7),end=currentdate),
window(stats::lag(urate,-5),start=c(1948,7),end=currentdate),
window(stats::lag(urate,-6),start=c(1948,7),end=currentdate)
)
colnames(uratelags)<-c("urate0","uratel1","uratel2","uratel3","uratel4",
"uratel5","uratel6")
ICSA<-fredr(series_id = "ICSA", vintage_dates=as.Date("2022-4-25")) #Initial Unemployment Claims, weekly
INDPRO <- fredr(series_id="INDPRO", units="chg", vintage_dates=as.Date("2022-4-25"))
# ICSA is weekly data: to make monthly, need to aggregate to lower frequency
xtsICSA<-xts(ICSA$value,order.by = ICSA$date) #Convert to irregular time series
#apply.monthly aggregates within a month using a function, here, the mean
claimsmonthly<-apply.monthly(xtsICSA,mean) #Can replace mean by other functions like max, var, etc
initialclaims<-ts(claimsmonthly,frequency=12,start=c(1967,1)) #convert back to monthly ts object
dinitialclaims<-diff(initialclaims) #Convert to month-over-month changes
indpro<-ts(INDPRO$value,frequency=12,start=c(1919,1),names="IndustrialProduction")
indprog<-window(diff(log(indpro)),start=c(1976,7))
Warning in log(indpro) : NaNs produced
xuratelags<-data.frame(
window(urate,start=c(2002,5),end=currentdate),
window(stats::lag(urate,-1),start=c(2002,5),end=currentdate),
window(stats::lag(urate,-2),start=c(2002,5),end=currentdate),
window(stats::lag(urate,-3),start=c(2002,5),end=currentdate),
window(stats::lag(urate,-4),start=c(2002,5),end=currentdate),
window(stats::lag(urate,-5),start=c(2002,5),end=currentdate),
window(stats::lag(urate,-6),start=c(2002,5),end=currentdate),
window(indprog,start=c(2002,5),end=currentdate),
window(dinitialclaims/1000,start=c(2002,5),end=currentdate)
)
colnames(xuratelags)<-c("urate0","uratel1","uratel2","uratel3","uratel4",
"uratel5","uratel6","indprog","dinitialclaims")
set.seed(12345)
utoday<-length(uratelags$urate0) #Last observation
ictoday<-length(dinitialclaims)
iptoday<-length(indprog)
xutodayslags<-data.frame(uratelags$urate0[utoday],
uratelags$uratel1[utoday],uratelags$uratel2[utoday],
uratelags$uratel3[utoday],uratelags$uratel4[utoday],
uratelags$uratel5[utoday],indprog[iptoday],dinitialclaims[ictoday]/1000)
names(xutodayslags)<-c("uratel1","uratel2","uratel3","uratel4","uratel5","uratel6","indprog","dinitialclaims")
maxorder<-6
defaultscale<-2.5
scalelist<-c(1,1,1,1,1,1,1,1) #6 lags, 2 external series
for (j in 1:maxorder){
scalelist[j]<-defaultscale*scalelist[j]/j^2
}
URATEbayesARX<-stan_glm(urate0~uratel1+uratel2+uratel3+uratel4+uratel5+
uratel6+indprog+dinitialclaims,
family=gaussian(),prior=normal(location=c(1,0,0,0,0,0,0.1,0.1),scale=scalelist),
data=xuratelags)
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 4.386 seconds (Warm-up)
Chain 1: 0.954 seconds (Sampling)
Chain 1: 5.34 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 3.799 seconds (Warm-up)
Chain 2: 0.992 seconds (Sampling)
Chain 2: 4.791 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.001 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 10 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 3.2 seconds (Warm-up)
Chain 3: 0.894 seconds (Sampling)
Chain 3: 4.094 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 3.531 seconds (Warm-up)
Chain 4: 0.999 seconds (Sampling)
Chain 4: 4.53 seconds (Total)
Chain 4:
#Sample from posterior predictive distribution
URATEbayesARXpp<-rstanarm::posterior_predict(URATEbayesARX,newdata=xutodayslags)
URATEbayesARXfcst<-mean(URATEbayesARXpp) #Forecast is posterior predictive mean
URATEbayesARXppint<-quantile(URATEbayesARXpp,c(0.025,0.975)) #Posterior predictive interval
plot(URATEbayesARX,"hist")
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ppdists<-data.frame(URATEbayesARXpp)
ggplot(data=ppdists)+geom_histogram(aes(x=X1))+
xlab("April Unemployment")+ylab("Posterior Frequency")+
ggtitle("Unemployment Posterior Predictions")
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
l1urate<-window(stats::lag(urate,-1),start=c(1949,1),end=c(2022,2))
l2urate<-window(stats::lag(urate,-2),start=c(1949,1),end=c(2022,2))
l3urate<-window(stats::lag(urate,-3),start=c(1949,1),end=c(2022,2))
l4urate<-window(stats::lag(urate,-4),start=c(1949,1),end=c(2022,2))
l5urate<-window(stats::lag(urate,-5),start=c(1949,1),end=c(2022,2))
l6urate<-window(stats::lag(urate,-6),start=c(1949,1),end=c(2022,2))
uraten<-window(urate,start=c(1949,1),end=c(2022,2))
udata<-data.frame(uraten,l1urate,l2urate,l3urate,l4urate,l5urate,l6urate)
udata12<-data.frame(uraten,l1urate)
options(mc.cores = parallel::detectCores())
mn <- bv_minnesota(lambda = bv_lambda(mode = 0.4), #tightness of prior covariance on coefficients
alpha = bv_alpha(mode = 2), #rate of decline of variance with increasing lag order
var = 1e06) #prior variance of constant term
mypriorsm <- bv_priors(hyper = c("lambda"), #Choose lambda by hierarchical method, leave others fixed
mn = mn) #Set Minnesota prior as defined above
min6u <- BVAR::bvar(udata,6,priors = mypriorsm)
Optimisation concluded.
Posterior marginal likelihood: 10452
Hyperparameters: lambda = 5
|
| | 0%
|
|= | 0%
|
|= | 1%
|
|== | 1%
|
|== | 2%
|
|=== | 2%
|
|=== | 3%
|
|==== | 3%
|
|==== | 4%
|
|===== | 4%
|
|===== | 5%
|
|====== | 5%
|
|====== | 6%
|
|======= | 6%
|
|======= | 7%
|
|======== | 7%
|
|======== | 8%
|
|========= | 8%
|
|========= | 9%
|
|========== | 9%
|
|========== | 10%
|
|=========== | 10%
|
|=========== | 11%
|
|============ | 11%
|
|============ | 12%
|
|============= | 12%
|
|============= | 13%
|
|============== | 13%
|
|============== | 14%
|
|=============== | 14%
|
|=============== | 15%
|
|================ | 15%
|
|================ | 16%
|
|================= | 16%
|
|================= | 17%
|
|================== | 17%
|
|================== | 18%
|
|=================== | 18%
|
|=================== | 19%
|
|==================== | 19%
|
|==================== | 20%
|
|===================== | 20%
|
|===================== | 21%
|
|====================== | 21%
|
|====================== | 22%
|
|======================= | 22%
|
|======================= | 23%
|
|======================== | 23%
|
|======================== | 24%
|
|========================= | 24%
|
|========================= | 25%
|
|========================== | 25%
|
|========================== | 26%
|
|=========================== | 26%
|
|=========================== | 27%
|
|============================ | 27%
|
|============================ | 28%
|
|============================= | 28%
|
|============================= | 29%
|
|============================== | 29%
|
|============================== | 30%
|
|=============================== | 30%
|
|=============================== | 31%
|
|================================ | 31%
|
|================================ | 32%
|
|================================= | 32%
|
|================================= | 33%
|
|================================== | 33%
|
|================================== | 34%
|
|=================================== | 34%
|
|=================================== | 35%
|
|==================================== | 35%
|
|==================================== | 36%
|
|===================================== | 36%
|
|===================================== | 37%
|
|====================================== | 37%
|
|====================================== | 38%
|
|======================================= | 38%
|
|======================================= | 39%
|
|======================================== | 39%
|
|======================================== | 40%
|
|========================================= | 40%
|
|========================================= | 41%
|
|========================================== | 41%
|
|========================================== | 42%
|
|=========================================== | 42%
|
|=========================================== | 43%
|
|============================================ | 43%
|
|============================================ | 44%
|
|============================================= | 44%
|
|============================================= | 45%
|
|============================================== | 45%
|
|============================================== | 46%
|
|=============================================== | 46%
|
|=============================================== | 47%
|
|================================================ | 47%
|
|================================================ | 48%
|
|================================================= | 48%
|
|================================================= | 49%
|
|================================================== | 49%
|
|================================================== | 50%
|
|=================================================== | 50%
|
|=================================================== | 51%
|
|==================================================== | 51%
|
|==================================================== | 52%
|
|===================================================== | 52%
|
|===================================================== | 53%
|
|====================================================== | 53%
|
|====================================================== | 54%
|
|======================================================= | 54%
|
|======================================================= | 55%
|
|======================================================== | 55%
|
|======================================================== | 56%
|
|========================================================= | 56%
|
|========================================================= | 57%
|
|========================================================== | 57%
|
|========================================================== | 58%
|
|=========================================================== | 58%
|
|=========================================================== | 59%
|
|============================================================ | 59%
|
|============================================================ | 60%
|
|============================================================= | 60%
|
|============================================================= | 61%
|
|============================================================== | 61%
|
|============================================================== | 62%
|
|=============================================================== | 62%
|
|=============================================================== | 63%
|
|================================================================ | 63%
|
|================================================================ | 64%
|
|================================================================= | 64%
|
|================================================================= | 65%
|
|================================================================== | 65%
|
|================================================================== | 66%
|
|=================================================================== | 66%
|
|=================================================================== | 67%
|
|==================================================================== | 67%
|
|==================================================================== | 68%
|
|===================================================================== | 68%
|
|===================================================================== | 69%
|
|====================================================================== | 69%
|
|====================================================================== | 70%
|
|======================================================================= | 70%
|
|======================================================================= | 71%
|
|======================================================================== | 71%
|
|======================================================================== | 72%
|
|========================================================================= | 72%
|
|========================================================================= | 73%
|
|========================================================================== | 73%
|
|========================================================================== | 74%
|
|=========================================================================== | 74%
|
|=========================================================================== | 75%
|
|============================================================================ | 75%
|
|============================================================================ | 76%
|
|============================================================================= | 76%
|
|============================================================================= | 77%
|
|============================================================================== | 77%
|
|============================================================================== | 78%
|
|=============================================================================== | 78%
|
|=============================================================================== | 79%
|
|================================================================================ | 79%
|
|================================================================================ | 80%
|
|================================================================================= | 80%
|
|================================================================================= | 81%
|
|================================================================================== | 81%
|
|================================================================================== | 82%
|
|=================================================================================== | 82%
|
|=================================================================================== | 83%
|
|==================================================================================== | 83%
|
|==================================================================================== | 84%
|
|===================================================================================== | 84%
|
|===================================================================================== | 85%
|
|====================================================================================== | 85%
|
|====================================================================================== | 86%
|
|======================================================================================= | 86%
|
|======================================================================================= | 87%
|
|======================================================================================== | 87%
|
|======================================================================================== | 88%
|
|========================================================================================= | 88%
|
|========================================================================================= | 89%
|
|========================================================================================== | 89%
|
|========================================================================================== | 90%
|
|=========================================================================================== | 90%
|
|=========================================================================================== | 91%
|
|============================================================================================ | 91%
|
|============================================================================================ | 92%
|
|============================================================================================= | 92%
|
|============================================================================================= | 93%
|
|============================================================================================== | 93%
|
|============================================================================================== | 94%
|
|=============================================================================================== | 94%
|
|=============================================================================================== | 95%
|
|================================================================================================ | 95%
|
|================================================================================================ | 96%
|
|================================================================================================= | 96%
|
|================================================================================================= | 97%
|
|================================================================================================== | 97%
|
|================================================================================================== | 98%
|
|=================================================================================================== | 98%
|
|=================================================================================================== | 99%
|
|==================================================================================================== | 99%
|
|==================================================================================================== | 100%
|
|=====================================================================================================| 100%
Finished MCMC after 20.96 secs.
bvar_preddu <- predict(min6u)
plot(bvar_preddu)
Y.sample1 <- window(uraten, end=c(2002, 12))
Y.sample2 <- window(uraten, start=c(2021,10))
end<-length(Y.sample2)
forecasts <- bvar_preddu[1]
forecastsw<-as.numeric(unlist(forecasts))[1:end]
rmse(Y.sample2, forecastsw)
[1] 0.367355
library(VARsignR)
Warning: package ‘VARsignR’ was built under R version 4.0.5
l1paymes2<-window(stats::lag(payems,-1),start=c(1949,1),end=c(2022,2))
l2paymes2<-window(stats::lag(payems,-2),start=c(1949,1),end=c(2022,2))
l3paymes2<-window(stats::lag(payems,-3),start=c(1949,1),end=c(2022,2))
l4paymes2<-window(stats::lag(payems,-4),start=c(1949,1),end=c(2022,2))
l5paymes2<-window(stats::lag(payems,-5),start=c(1949,1),end=c(2022,2))
l6paymes2<-window(stats::lag(payems,-6),start=c(1949,1),end=c(2022,2))
payemsn2<-window(payems,start=c(1949,1),end=c(2022,2))
pdata2<-data.frame(payemsn2,l1paymes2,l2paymes2,l3paymes2,l4paymes2,l5paymes2,l6paymes2)
jdat<-data.frame(udata,pdata2)
jdat<-ts(jdat)
jdata12<-data.frame(uraten,payemsn2)
jdata12<-ts(jdata12)
library(bvartools)
Warning: package ‘bvartools’ was built under R version 4.0.5
Loading required package: coda
Warning: package ‘coda’ was built under R version 4.0.5
Attaching package: ‘coda’
The following object is masked from ‘package:Boom’:
thin
Registered S3 methods overwritten by 'bvartools':
method from
plot.bvar BVAR
predict.bvar BVAR
summary.bvar BVAR
Attaching package: ‘bvartools’
The following objects are masked from ‘package:BVAR’:
bvar, fevd, irf
The following objects are masked from ‘package:vars’:
fevd, irf
nwu<- gen_var(jdata12, p = 2, deterministic = "const",
iterations = 5000, burnin = 1000)
nwup <-add_priors(nwu,
coef = list(v_i = 1, v_i_det = 0.1, shape = 3, rate = 1e-04, rate_det = 0.01),
coint = list(v_i = 0, p_tau_i = 1, shape = 3, rate = 1e-04, rho = 0.999),
sigma = list(df = "k", scale = 1, mu = 0, v_i = 0.01, sigma_h = 0.05),
ssvs = NULL,
bvs = NULL)
bvar_est <- draw_posterior(nwup)
Estimating model...
bvar_pred <- predict(bvar_est, n.ahead = 10, new_d = rep(1, 10))
plot(bvar_pred)
##dirichlet prior bvar omitted
ICSA<-fredr(series_id = "ICSA", vintage_dates=as.Date("2022-3-22")) #Initial Unemployment Claims, weekly
NPPTTL<-fredr(series_id = "NPPTTL",units="chg",vintage_dates=as.Date("2022-3-22")) #ADP Payrolls survey
# ICSA is weekly data: to make monthly, need to aggregate to lower frequency
xtsICSA<-xts(ICSA$value,order.by = ICSA$date) #Convert to irregular time series
#apply.monthly aggregates within a month using a function, here, the mean
claimsmonthly<-apply.monthly(xtsICSA,mean) #Can replace mean by other functions like max, var, etc
initialclaims<-ts(claimsmonthly,frequency=12,start=c(1967,1)) #convert back to monthly ts object
dinitialclaims<-diff(initialclaims) #Convert to month-over-month changes
#Convert to ts object, get rid of initial NA values
nppttl<-window(ts(NPPTTL$value,frequency=12,start=c(2000,12)),start=c(2002,5))
payemslags<-data.frame(
window(payems,start=c(1939,8),end=currentdate),
window(stats::lag(payems,-1),start=c(1939,8),end=currentdate),
window(stats::lag(payems,-2),start=c(1939,8),end=currentdate),
window(stats::lag(payems,-3),start=c(1939,8),end=currentdate),
window(stats::lag(payems,-4),start=c(1939,8),end=currentdate),
window(stats::lag(payems,-5),start=c(1939,8),end=currentdate),
window(stats::lag(payems,-6),start=c(1939,8),end=currentdate)
)
colnames(payemslags)<-c("payems0","payemsl1","payemsl2","payemsl3","payemsl4","payemsl5","payemsl6")
xuratelags<-data.frame(
window(urate,start=c(2002,6),end=currentdate),
window(stats::lag(urate,-1),start=c(2002,6),end=currentdate),
window(stats::lag(urate,-2),start=c(2002,6),end=currentdate),
window(stats::lag(urate,-3),start=c(2002,6),end=currentdate),
window(stats::lag(urate,-4),start=c(2002,6),end=currentdate),
window(stats::lag(urate,-5),start=c(2002,6),end=currentdate),
window(stats::lag(urate,-6),start=c(2002,6),end=currentdate),
window(nppttl/100,start=c(2002,5),end=currentdate),
window(dinitialclaims/1000,start=c(2002,6),end=currentdate)
)
Warning in window.default(x, ...) : 'end' value not changed
colnames(xuratelags)<-c("urate0","uratel1","uratel2","uratel3","uratel4",
"uratel5","uratel6","nppttl","dinitialclaims")
xpayemslags<-data.frame(
window(payems,start=c(2002,6),end=currentdate),
window(stats::lag(payems,-1),start=c(2002,6),end=currentdate),
window(stats::lag(payems,-2),start=c(2002,6),end=currentdate),
window(stats::lag(payems,-3),start=c(2002,6),end=currentdate),
window(stats::lag(payems,-4),start=c(2002,6),end=currentdate),
window(stats::lag(payems,-5),start=c(2002,6),end=currentdate),
window(stats::lag(payems,-6),start=c(2002,6),end=currentdate),
window(nppttl,start=c(2002,5),end=currentdate),
window(dinitialclaims/1000,start=c(2002,6),end=currentdate)
)
Warning in window.default(x, ...) : 'end' value not changed
colnames(xpayemslags)<-c("payems0","payemsl1","payemsl2","payemsl3",
"payemsl4","payemsl5","payemsl6","nppttl","dinitialclaims")
#Construct current month info to build predictions
utoday<-length(uratelags$urate0) #Last observation
adptoday<-length(nppttl)
ictoday<-length(dinitialclaims)
xutodayslags<-data.frame(uratelags$urate0[utoday],
uratelags$uratel1[utoday],uratelags$uratel2[utoday],
uratelags$uratel3[utoday],uratelags$uratel4[utoday],
uratelags$uratel5[utoday],nppttl[adptoday]/100,dinitialclaims[ictoday]/1000)
ptoday<-length(payemslags$payems0) #Last observation
xptodayslags<-data.frame(payemslags$payems0[ptoday],payemslags$payemsl1[ptoday],
payemslags$payemsl2[ptoday],payemslags$payemsl3[ptoday],payemslags$payemsl4[ptoday],
payemslags$payemsl5[ptoday],nppttl[adptoday],dinitialclaims[ictoday]/1000)
#Rename everything
names(xutodayslags)<-c("uratel1","uratel2","uratel3","uratel4","uratel5","uratel6","nppttl","dinitialclaims")
names(xptodayslags)<-c("payemsl1","payemsl2","payemsl3","payemsl4","payemsl5","payemsl6","nppttl","dinitialclaims")
xmatrix<-as.matrix(data.frame( window(dinitialclaims,start=c(2002,6),end=currentdate)/1000,
window(nppttl,start=c(2002,5),end=currentdate)) )
Warning in window.default(x, ...) : 'end' value not changed
colnames(xmatrix)<-c("dinitclaims","nppttl")
uratec<-window(urate,start=c(2002,6))
payemsc<-window(payems,start=c(2002,6))
URATEARIMAX<-auto.arima(uratec,xreg=xmatrix,seasonal=FALSE)
PAYEMSARIMAX<-auto.arima(payemsc,xreg=xmatrix,seasonal=FALSE)
initctoday<-dinitialclaims[ictoday]/1000
nppttloday<-nppttl[adptoday]
xtoday=as.matrix(data.frame(initctoday,nppttloday))
colnames(xtoday)<-c("dinitclaims","nppttl")
URATEARIMAXfcst <-forecast::forecast(URATEARIMAX,xreg=xtoday,h=1)
PAYEMSARIMAXfcst <-forecast::forecast(PAYEMSARIMAX,xreg=xtoday,h=1)
PAYEMSARIMAXfcst
PRevious and Baseline Forecasts
options(mc.cores = parallel::detectCores()) #speed up sampling
set.seed(123546) #make results reproducible
bayesuratec<-stan_glm(uraten~l1urate+l2urate+l3urate+l4urate+l5urate+l6urate,
data=udata,
family = gaussian(),
prior=laplace(location=c(1,0,0,0,0,0),scale=c(2.5,2.5/(2^2),2.5/(3^2),2.5/(4^2),2.5/(5^2),2.5/(6^2))),
prior_intercept=normal(location=0, scale=10)
)
l1paymes<-window(stats::lag(payems,-1),start=c(1940,1),end=c(2022,2))
l2paymes<-window(stats::lag(payems,-2),start=c(1940,1),end=c(2022,2))
l3paymes<-window(stats::lag(payems,-3),start=c(1940,1),end=c(2022,2))
l4paymes<-window(stats::lag(payems,-4),start=c(1940,1),end=c(2022,2))
l5paymes<-window(stats::lag(payems,-5),start=c(1940,1),end=c(2022,2))
l6paymes<-window(stats::lag(payems,-6),start=c(1940,1),end=c(2022,2))
payemsn<-window(payems,start=c(1940,1),end=c(2022,2))
pdata<-data.frame(payemsn,l1paymes,l2paymes,l3paymes,l4paymes,l5paymes,l6paymes)
min6p <- BVAR::bvar(pdata,6,priors = mypriorsm)
Optimisation concluded.
Posterior marginal likelihood: -6611.11
Hyperparameters: lambda = 5
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|=== | 3%
|
|=== | 4%
|
|==== | 4%
|
|==== | 5%
|
|===== | 5%
|
|===== | 6%
|
|====== | 6%
|
|====== | 7%
|
|======= | 7%
|
|======= | 8%
|
|======== | 8%
|
|======== | 9%
|
|========= | 9%
|
|========= | 10%
|
|========== | 10%
|
|========== | 11%
|
|=========== | 11%
|
|=========== | 12%
|
|============ | 12%
|
|============ | 13%
|
|============ | 14%
|
|============= | 14%
|
|============= | 15%
|
|============== | 15%
|
|============== | 16%
|
|=============== | 16%
|
|=============== | 17%
|
|================ | 17%
|
|================ | 18%
|
|================= | 18%
|
|================= | 19%
|
|================== | 19%
|
|================== | 20%
|
|=================== | 20%
|
|=================== | 21%
|
|==================== | 21%
|
|==================== | 22%
|
|===================== | 22%
|
|===================== | 23%
|
|====================== | 23%
|
|====================== | 24%
|
|======================= | 24%
|
|======================= | 25%
|
|======================= | 26%
|
|======================== | 26%
|
|======================== | 27%
|
|========================= | 27%
|
|========================= | 28%
|
|========================== | 28%
|
|========================== | 29%
|
|=========================== | 29%
|
|=========================== | 30%
|
|============================ | 30%
|
|============================ | 31%
|
|============================= | 31%
|
|============================= | 32%
|
|============================== | 32%
|
|============================== | 33%
|
|=============================== | 33%
|
|=============================== | 34%
|
|================================ | 34%
|
|================================ | 35%
|
|================================= | 35%
|
|================================= | 36%
|
|================================== | 36%
|
|================================== | 37%
|
|================================== | 38%
|
|=================================== | 38%
|
|=================================== | 39%
|
|==================================== | 39%
|
|==================================== | 40%
|
|===================================== | 40%
|
|===================================== | 41%
|
|====================================== | 41%
|
|====================================== | 42%
|
|======================================= | 42%
|
|======================================= | 43%
|
|======================================== | 43%
|
|======================================== | 44%
|
|========================================= | 44%
|
|========================================= | 45%
|
|========================================== | 45%
|
|========================================== | 46%
|
|=========================================== | 46%
|
|=========================================== | 47%
|
|============================================ | 47%
|
|============================================ | 48%
|
|============================================= | 48%
|
|============================================= | 49%
|
|============================================== | 49%
|
|============================================== | 50%
|
|============================================== | 51%
|
|=============================================== | 51%
|
|=============================================== | 52%
|
|================================================ | 52%
|
|================================================ | 53%
|
|================================================= | 53%
|
|================================================= | 54%
|
|================================================== | 54%
|
|================================================== | 55%
|
|=================================================== | 55%
|
|=================================================== | 56%
|
|==================================================== | 56%
|
|==================================================== | 57%
|
|===================================================== | 57%
|
|===================================================== | 58%
|
|====================================================== | 58%
|
|====================================================== | 59%
|
|======================================================= | 59%
|
|======================================================= | 60%
|
|======================================================== | 60%
|
|======================================================== | 61%
|
|========================================================= | 61%
|
|========================================================= | 62%
|
|========================================================== | 62%
|
|========================================================== | 63%
|
|========================================================== | 64%
|
|=========================================================== | 64%
|
|=========================================================== | 65%
|
|============================================================ | 65%
|
|============================================================ | 66%
|
|============================================================= | 66%
|
|============================================================= | 67%
|
|============================================================== | 67%
|
|============================================================== | 68%
|
|=============================================================== | 68%
|
|=============================================================== | 69%
|
|================================================================ | 69%
|
|================================================================ | 70%
|
|================================================================= | 70%
|
|================================================================= | 71%
|
|================================================================== | 71%
|
|================================================================== | 72%
|
|=================================================================== | 72%
|
|=================================================================== | 73%
|
|==================================================================== | 73%
|
|==================================================================== | 74%
|
|===================================================================== | 74%
|
|===================================================================== | 75%
|
|===================================================================== | 76%
|
|====================================================================== | 76%
|
|====================================================================== | 77%
|
|======================================================================= | 77%
|
|======================================================================= | 78%
|
|======================================================================== | 78%
|
|======================================================================== | 79%
|
|========================================================================= | 79%
|
|========================================================================= | 80%
|
|========================================================================== | 80%
|
|========================================================================== | 81%
|
|=========================================================================== | 81%
|
|=========================================================================== | 82%
|
|============================================================================ | 82%
|
|============================================================================ | 83%
|
|============================================================================= | 83%
|
|============================================================================= | 84%
|
|============================================================================== | 84%
|
|============================================================================== | 85%
|
|=============================================================================== | 85%
|
|=============================================================================== | 86%
|
|================================================================================ | 86%
|
|================================================================================ | 87%
|
|================================================================================ | 88%
|
|================================================================================= | 88%
|
|================================================================================= | 89%
|
|================================================================================== | 89%
|
|================================================================================== | 90%
|
|=================================================================================== | 90%
|
|=================================================================================== | 91%
|
|==================================================================================== | 91%
|
|==================================================================================== | 92%
|
|===================================================================================== | 92%
|
|===================================================================================== | 93%
|
|====================================================================================== | 93%
|
|====================================================================================== | 94%
|
|======================================================================================= | 94%
|
|======================================================================================= | 95%
|
|======================================================================================== | 95%
|
|======================================================================================== | 96%
|
|========================================================================================= | 96%
|
|========================================================================================= | 97%
|
|========================================================================================== | 97%
|
|========================================================================================== | 98%
|
|=========================================================================================== | 98%
|
|=========================================================================================== | 99%
|
|============================================================================================| 99%
|
|============================================================================================| 100%
Finished MCMC after 19.86 secs.
#bvar_preddp<- stats::predict(min6p)
#plot(bvar_preddp)
INDPRO <- fredr(series_id = "INDPRO")
tIND <- ts(INDPRO$value,frequency=12,start=c(1959,1),end = c(2022,1),names="Industrial Production")
M1<-fredr(series_id = "M1SL")
m1<-ts(M1$value,frequency = 12,start=c(1959,1),names="M1 Growth")
jointdata2 <- ts(data.frame("urate" = window(urate, start = c(1959,3)),"payems" = window(payems, start = c(1959,3))))
varinc1<-VAR(jointdata2, p = 1)
varinc2<-VAR(jointdata2, p = 2)
varinc3<-VAR(jointdata2, p = 3)
varinc4<-VAR(jointdata2, p = 4)
varinc5<-VAR(jointdata2, p = 5)
jdvar1 <-forecast(varinc1, h = 4)
jdvar2 <-forecast(varinc2, h = 4)
jdvar3 <-forecast(varinc3, h = 4)
jdvar4 <-forecast(varinc4, h = 4)
jdvar5 <-forecast(varinc5, h = 4)
jdvar1
urate
payems
jdvar2
urate
payems
jdvar3
urate
payems
jdvar4
urate
payems
jdvar5
urate
payems
jdvar<-data.frame(jdvar1,jdvar2,jdvar3,jdvar4,jdvar5)
autoplot(jdvar1)
autoplot(jdvar2)
autoplot(jdvar3)
autoplot(jdvar4)
autoplot(jdvar5)
naive(payems)
unemAR1selection<-Arima(urate,order=c(1,0,0),include.mean=FALSE)
payemsAR1selection<-Arima(payems,order=c(1,0,0),include.mean=FALSE)
unemAR1<-forecast(unemAR1selection,h=1)
payemsAR1<-forecast(payemsAR1selection,h=1)
unemAR1
payemsAR1
plot(unemAR1)
plot(payemsAR1)
far1 <- function(x, h){forecast(Arima(x, order=c(1,0,0)), h=h)}
far2 <- function(x, h){forecast(Arima(x, order=c(2,0,0)), h=h)}
far3 <- function(x, h){forecast(Arima(x, order=c(3,0,0)), h=h)}
far4 <- function(x, h){forecast(Arima(x, order=c(4,0,0)), h=h)}
far5 <- function(x, h){forecast(Arima(x, order=c(5,0,0)), h=h)}
urateAR1<- far1(urate, 4)
urateAR2<- far2(urate, 4)
urateAR3<- far3(urate, 4)
urateAR4<- far4(urate, 4)
urateAR5<- far5(urate, 4)
urateAR1
urateAR2
urateAR3
urateAR4
urateAR5
autoplot(urateAR1)+
ggtitle("AR1 forecast for unemployment") +
xlab("Year") +
ylab("Unemployment Rate")
autoplot(urateAR2)+
ggtitle("AR2 forecast for unemployment") +
xlab("Year") +
ylab("Unemployment Rate")
autoplot(urateAR3)+
ggtitle("AR3 forecast for unemployment") +
xlab("Year") +
ylab("Unemployment Rate")
autoplot(urateAR4)+
ggtitle("AR4 forecast for unemployment") +
xlab("Year") +
ylab("Unemployment Rate")
autoplot(urateAR5)+
ggtitle("AR5 forecast for unemployment") +
xlab("Year") +
ylab("Unemployment Rate")
payemsAR1 <- far1(payems, 4)
payemsAR2 <- far2(payems, 4)
payemsAR3 <- far3(payems, 4)
payemsAR4 <- far4(payems, 4)
payemsAR5 <- far4(payems, 4)
payemsAR1
payemsAR2
payemsAR3
payemsAR4
payemsAR5
autoplot(payemsAR1)+
ggtitle("AR1 forecast fornon farm payroll growth") +
xlab("Year") +
ylab("Payroll Growth")
Warning: Removed 1 row(s) containing missing values (geom_path).
autoplot(payemsAR2)+
ggtitle("AR2 forecast fornon farm payroll growth") +
xlab("Year") +
ylab("Payroll Growth")
Warning: Removed 1 row(s) containing missing values (geom_path).
autoplot(payemsAR3)+
ggtitle("AR3 forecast fornon farm payroll growth") +
xlab("Year") +
ylab("Payroll Growth")
Warning: Removed 1 row(s) containing missing values (geom_path).
autoplot(payemsAR4)+
ggtitle("AR4 forecast fornon farm payroll growth") +
xlab("Year") +
ylab("Payroll Growth")
Warning: Removed 1 row(s) containing missing values (geom_path).
autoplot(payemsAR5)+
ggtitle("AR5 forecast fornon farm payroll growth") +
xlab("Year") +
ylab("Payroll Growth")
Warning: Removed 1 row(s) containing missing values (geom_path).
uratel1 <- window(diff(urate, lag = 1), start = c(1948,6))
uratel2 <- window(diff(urate, lag = 2), start = c(1948,6))
uratel3 <- window(diff(urate, lag = 3), start = c(1948,6))
uratel4 <- window(diff(urate, lag = 4), start = c(1948,6))
uratel5 <- window(diff(urate, lag = 5), start = c(1948,6))
uraten <- window(urate, start = c(1948,6))
Unemlag <- data.frame(cbind(uraten,uratel1,uratel2,uratel3,uratel4,uratel5))
rqu1 <- rq(uraten~uratel1)
rqu2 <- rq(uraten~uratel1 + uratel2)
rqu3 <- rq(uraten~uratel1 + uratel2 + uratel3)
rqu4 <- rq(uraten~uratel1 + uratel2 + uratel3 + uratel4)
rqu5 <- rq(uraten~uratel1 + uratel2 + uratel3 + uratel4 + uratel5)
uARf1<-predict(rqu1,Unemlag)
uARf2<-predict(rqu2,Unemlag)
uARf3<-predict(rqu3,Unemlag)
uARf4<-predict(rqu4,Unemlag)
uARf5<-predict(rqu5,Unemlag)
payemsl1 <- window(diff(payems, lag = 1), start = c(1948,6))
payemsl2 <- window(diff(payems, lag = 2), start = c(1948,6))
payemsl3 <- window(diff(payems, lag = 3), start = c(1948,6))
payemsl4 <- window(diff(payems, lag = 4), start = c(1948,6))
payemsl5 <- window(diff(payems, lag = 5), start = c(1948,6))
payemsn <- window(payems, start = c(1948,6))
Payemslag <- data.frame(cbind(payemsn,payemsl1,payemsl2,payemsl3,payemsl4,payemsl5))
rqp1 <- rq(payemsn~payemsl1)
rqp2 <- rq(payemsn~payemsl1 + payemsl2)
rqp3 <- rq(payemsn~payemsl1 + payemsl2 + payemsl3)
rqp4 <- rq(payemsn~payemsl1 + payemsl2 + payemsl3 + payemsl4)
rqp5 <- rq(payemsn~payemsl1 + payemsl2 + payemsl3 + payemsl4 + payemsl5)
pARf1<- predict(rqp1,Payemslag)
pARf2<-predict(rqp2,Payemslag)
pARf3<-predict(rqp3,Payemslag)
pARf4<-predict(rqp4,Payemslag)
pARf5<-predict(rqp5,Payemslag)
forecast::accuracy(urateAR1)[2]
[1] 0.4177033
forecast::accuracy(urateAR2)[2]
[1] 0.4170023
forecast::accuracy(urateAR3)[2]
[1] 0.4169268
forecast::accuracy(urateAR4)[2]
[1] 0.4166857
forecast::accuracy(urateAR5)[2]
[1] 0.4163575
forecast::accuracy(payemsAR1)[3]
[1] 192.002
forecast::accuracy(payemsAR2)[3]
[1] 202.2887
forecast::accuracy(payemsAR3)[3]
[1] 201.6351
forecast::accuracy(payemsAR4)[3]
[1] 203.7887
forecast::accuracy(payemsAR5)[3]
[1] 203.7887
mean(abs(uraten - uARf1))
[1] 1.293266
mean(abs(uraten - uARf2))
[1] 1.286682
mean(abs(uraten - uARf3))
[1] 1.278053
mean(abs(uraten - uARf4))
[1] 1.269127
mean(abs(uraten - uARf5))
[1] 1.259105
mean(abs(payemsn - pARf1))
[1] 169.5133
mean(abs(payemsn - pARf2))
[1] 158.3281
mean(abs(payemsn - pARf3))
[1] 150.873
mean(abs(payemsn - pARf4))
[1] 145.3968
mean(abs(payemsn - pARf5))
[1] 140.8658
forecast::accuracy(jdvar1, d=0,D = 0)
ME RMSE MAE MPE MAPE MASE
urate Training set -5.928105e-17 0.4311465 0.1630065 -0.3831136 2.637936 0.02729032
payems Training set 5.872287e-15 802.2704926 200.8587910 73.2644574 182.163661 0.79546996
ACF1
urate Training set -0.004309431
payems Training set 0.006962061
forecast::accuracy(jdvar2, d=0,D = 0)
ME RMSE MAE MPE MAPE MASE
urate Training set 2.886950e-17 0.4295283 0.1628046 -0.3726559 2.632027 0.02725652
payems Training set -4.730579e-16 795.0658943 208.0935341 88.3776064 202.000374 0.82412203
ACF1
urate Training set 0.002381700
payems Training set -0.003404134
forecast::accuracy(jdvar3, d=0,D = 0)
ME RMSE MAE MPE MAPE MASE
urate Training set -2.376761e-18 0.4289234 0.1644295 -0.3722356 2.661296 0.02752856
payems Training set 1.596755e-15 794.0064990 210.8322539 85.1072029 211.490746 0.83496830
ACF1
urate Training set -0.000939372
payems Training set -0.002200142
forecast::accuracy(jdvar4, d=0,D = 0)
ME RMSE MAE MPE MAPE MASE
urate Training set -5.956309e-17 0.4270461 0.1658112 -0.3625076 2.686989 0.02775989
payems Training set 1.698630e-14 793.4365962 212.0268720 89.1909360 216.103299 0.83969940
ACF1
urate Training set -0.0015806782
payems Training set -0.0005104153
forecast::accuracy(jdvar5, d=0,D = 0)
ME RMSE MAE MPE MAPE MASE
urate Training set -1.006239e-17 0.4267696 0.1657163 -0.3614512 2.688204 0.02774399
payems Training set -4.097181e-14 793.5520448 211.1353430 88.6947865 216.897871 0.83616864
ACF1
urate Training set -0.002156346
payems Training set -0.000950823
uratenw<-window(uraten,start=c(1949,1),end=c(2022,2))
historyu <- data.frame(ds = seq(as.Date('1949-01-01'), as.Date('2022-02-01'), by = 'm'), y = uratenw)
myprophetmodelu<-prophet(
df = historyu,
growth = "linear",
changepoints = NULL,
n.changepoints = 25,
changepoint.range = 0.8,
yearly.seasonality = "auto",
weekly.seasonality = "auto",
daily.seasonality = "auto",
holidays = NULL,
seasonality.mode = "additive",
seasonality.prior.scale = 10,
holidays.prior.scale = 10,
changepoint.prior.scale = 0.05,
mcmc.samples = 0,
interval.width = 0.8,
uncertainty.samples = 1000,
fit = TRUE
)
Disabling weekly seasonality. Run prophet with weekly.seasonality=TRUE to override this.
Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
#Produce forecasts at monthly frequency
futureu <- make_future_dataframe(myprophetmodelu, periods = 4,freq = "month")
myprophetfcstu<-predict(myprophetmodelu,futureu)
#Command to plot forecast (uses both fit object and forecast object)
plot(myprophetmodelu,myprophetfcstu)
payemsnw<-window(payemsn,end=c(2022,2))
historyp <- data.frame(ds = seq(as.Date('1948-06-01'), as.Date('2022-02-01'), by = 'm'), y = payemsnw)
myprophetmodelp<-prophet(
df = historyp,
growth = "linear",
changepoints = NULL,
n.changepoints = 25,
changepoint.range = 0.8,
yearly.seasonality = "auto",
weekly.seasonality = "auto",
daily.seasonality = "auto",
holidays = NULL,
seasonality.mode = "additive",
seasonality.prior.scale = 10,
holidays.prior.scale = 10,
changepoint.prior.scale = 0.05,
mcmc.samples = 0,
interval.width = 0.8,
uncertainty.samples = 1000,
fit = TRUE
)
Disabling weekly seasonality. Run prophet with weekly.seasonality=TRUE to override this.
Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
#Produce forecasts at monthly frequency
futurep <- make_future_dataframe(myprophetmodelp, periods = 4,freq = "month")
myprophetfcstp<-predict(myprophetmodelp,futurep)
#Command to plot forecast (uses both fit object and forecast object)
plot(myprophetmodelp,myprophetfcstp)