1 Introduction

This quick note is intended to introduce the intuition behind the NNS.VAR function, which serves as a nonparametric vector autoregression.

Install latest version of NNS (>=0.4.7) and other required packages…

require(devtools)
install_github("OVVO-Financial/NNS", ref = "NNS-Beta-Version")

library(NNS)
library(vars)
library(forecast)
library(randomForest)

library(kableExtra)

2 `NNS.VAR` Example with 3 Variables: Fed Funds, Real GNP, and Inflation

Using the log differences of each Real GNP and Inflation, the previous method creates lowest frequency series. In this illustrative VAR example from Lutz Kilian¹, Fed Funds (monthly) is averaged over 3-period windows to align with quarterly GNP data.

# Load Variables
FF = read.csv("fedfunds.txt", header = FALSE, sep = "\t")

realgnp = read.csv("realgnp.txt", header = FALSE, sep = "\t")
drgdp = diff(log(realgnp[, 3])) * 100

gnpdeflator = read.csv("gnpdeflator.txt", header = FALSE, sep = "\t")
infl = diff(log(gnpdeflator[, 3])) * 100

irate = numeric()
for (i in seq(1, length(FF[, 3]), 3)) {
    irate[i] = mean(FF[i:(i + 2), 3])
}

irate = na.omit(irate)


y = cbind(`Diff Real GDP` = drgdp[1:213], `Interest Rate` = irate, 
    Inflation = infl[1:213])

head(y)

##      Diff Real GDP Interest Rate Inflation
## [1,]     1.9944940     0.9866667 0.2810368
## [2,]     2.8100061     1.3433333 0.4818031
## [3,]     1.5928111     1.5000000 0.4148308
## [4,]     1.3345828     1.9400000 0.7089943
## [5,]     0.6089062     2.3566667 0.9778623
## [6,]    -0.3111535     2.4833333 0.9998818

tail(y)

##        Diff Real GDP Interest Rate Inflation
## [208,]   -0.01158863      5.246667 0.6891565
## [209,]    0.84719231      5.246667 0.3501309
## [210,]    0.06757442      5.256667 1.1124581
## [211,]    0.91316852      5.250000 0.5562536
## [212,]    1.05338427      5.073333 0.3451964
## [213,]    0.66741511      4.496667 0.4471854

3 Forecast 12 Quarters Out-of-Sample Using All Methods

We will withhold 12 quarters of observations as the test set. Let’s see how accurate a VAR, a random forest and NNS are…

Why a random forest? The Bank of England has recently been posting working papers utilizing machine learning, specifically random forests:

4 Traditional VAR

# Forecast next 12 quarters
h = 12

# Create train and test sets for both versions
y_train_VAR = head(y, dim(y)[1] - h)
y_test_VAR = tail(y, h)

# 4 quarterly lags used in VAR
VAR_train = VAR(y = y_train_VAR, p = 4)
VAR_test = predict(VAR_train, n.ahead = h)

par(mfrow = c(1, 3))


VAR_predictions_RMSE = cbind(rbind(Diff.Real.GDP = sqrt(mean((VAR_test$fcst[[1]][, 
    1] - y_test_VAR[, 1])^2)), Interest.Rate = sqrt(mean((VAR_test$fcst[[2]][, 
    1] - y_test_VAR[, 2])^2)), Inflation = sqrt(mean((VAR_test$fcst[[3]][, 
    1] - y_test_VAR[, 3])^2))))

colnames(VAR_predictions_RMSE) = c("VAR Estimates RMSE to Actual Values")


VAR_predictions_corr = cbind(rbind(Diff.Real.GDP = cor(VAR_test$fcst[[1]][, 
    1], y_test_VAR[, 1], method = "spearman"), Interest.Rate = cor(VAR_test$fcst[[2]][, 
    1], y_test_VAR[, 2], method = "spearman"), Inflation = cor(VAR_test$fcst[[3]][, 
    1], y_test_VAR[, 3], method = "spearman")))

colnames(VAR_predictions_corr) = c("VAR Estimates Correlation to Actual Values")

VAR_results = cbind(VAR_predictions_RMSE, VAR_predictions_corr)
knitr::kable(VAR_results, digits = 4) %>% kable_styling(full_width = T) %>% 
    column_spec(1, width = "8cm")

	VAR Estimates RMSE to Actual Values	VAR Estimates Correlation to Actual Values
Diff.Real.GDP	0.5062	-0.3916
Interest.Rate	1.3451	0.7741
Inflation	0.3481	-0.6294

The correlations for the VAR forecast of realgdp and inflation versus the actual out-of-sample values were negative! The ARMA forecasts were negative for realgdp and interest rates.

This highlights the dual objective of differenced time series data:

Get the sign right
Get the magnitude right

NNS is able to address this dual objective into its routines as detailed in the following code chunks.

5 `NNS` and Random Forest

5.1 Function to Generate Lagged Variables

The following code is to generate a matrix of lagged variables to be used in both NNS and Random Forest routines.

lag.mtx <- function(x, tau) {
    colheads <- NULL
    
    if (is.null(dim(x)[2])) {
        colheads <- noquote(as.character(deparse(substitute(x))))
        x <- t(t(x))
    }
    
    j.vectors <- list()
    
    for (j in 1:ncol(x)) {
        if (is.null(colheads)) {
            colheads <- colnames(x)[j]
            
            colheads <- noquote(as.character(deparse(substitute(colheads))))
        }
        
        x.vectors <- list()
        heads <- paste0(colheads, ".tau.")
        heads <- gsub("\"", "", heads)
        
        for (i in 0:tau) {
            x.vectors[[paste(heads, i, sep = "")]] <- numeric(0L)
            start <- tau - i + 1
            end <- length(x[, j]) - i
            x.vectors[[i + 1]] <- x[start:end, j]
        }
        
        j.vectors[[j]] <- do.call(cbind, x.vectors)
        colheads <- NULL
    }
    
    return(as.data.frame(do.call(cbind, j.vectors)))
}


# Test it
set.seed(123)
V1 <- rnorm(10)
V2 <- rnorm(10)
V3 <- rnorm(10)

lag.mtx(cbind(V1 = V1, V2 = V2, V3 = V3), 3)

# Compare to last values of each variable
cbind(V1 = tail(V1), V2 = tail(V2), V3 = tail(V3))

##              V1         V2         V3
## [1,]  0.1292877 -0.5558411 -0.6250393
## [2,]  1.7150650  1.7869131 -1.6866933
## [3,]  0.4609162  0.4978505  0.8377870
## [4,] -1.2650612 -1.9666172  0.1533731
## [5,] -0.6868529  0.7013559 -1.1381369
## [6,] -0.4456620 -0.4727914  1.2538149

6 Predictions

6.1 1 Year Lag \((\tau = 4)\), 3 Year Forecast \((h = 12)\)

6.2 Create Forecasted IVs Using `NNS.ARMA()`

train_l = dim(y)[1] - h
test_DVs = tail(y, h)

nns_IVs = list()

cl <- makeCluster(detectCores() - 1)
registerDoParallel(cl)

nns_IVs <- foreach(i = 1:ncol(y_train_VAR), .packages = "NNS") %dopar% 
    {
        variable = y_train_VAR[, i]
        
        periods = NNS.seas(variable, modulo = 4, mod.only = FALSE, 
            plot = FALSE)$periods
        
        b = NNS.ARMA.optim(variable, seasonal.factor = periods, 
            training.set = length(variable) - 2 * h, ncores = 1, 
            print.trace = FALSE, obj.fn = expression(sum((predicted - 
                actual)^2)), objective = "min")
        
        NNS.ARMA(variable, h = h, seasonal.factor = b$periods, 
            weights = b$weights, method = b$method, ncores = 1, 
            plot = FALSE) + b$bias.shift
    }

stopCluster(cl)
registerDoSEQ()

nns_IVs = do.call(cbind, nns_IVs)


# Combine forecasted IVs onto training data.frame
new_values = rbind(y_train_VAR, nns_IVs)

# Now lag new forecasted data.frame
lagged_new_values = lag.mtx(new_values, tau = 4)


# Test how accurate new univariate NNS.ARMA forecast IVs are
# Also add auto arima estimate...
NNS.ARMA_RMSEs = numeric()
NNS.ARMA_Correlations = numeric()
auto.ARMA_RMSEs = numeric()
auto.ARMA_Correlations = numeric()

arma.fit = list()

for (i in 1:3) {
    NNS.ARMA_RMSEs[i] = sqrt(mean((nns_IVs[, i] - test_DVs[, 
        i])^2))
    NNS.ARMA_Correlations[i] = cor(nns_IVs[, i], test_DVs[, i], 
        method = "spearman")
    fit = auto.arima(y_train_VAR[, i])
    arma.fit[[i]] = as.numeric(forecast(fit, h = h)$mean)
    auto.ARMA_RMSEs[i] = sqrt(mean((arma.fit[[i]] - test_DVs[, 
        i])^2))
    auto.ARMA_Correlations[i] = cor(arma.fit[[i]], test_DVs[, 
        i], method = "spearman")
}

NNS.ARMA_results = cbind(NNS.ARMA_RMSEs, NNS.ARMA_Correlations)

colnames(NNS.ARMA_results) = c("NNS ARMA RMSEs", "NNS ARMA Corr")

knitr::kable(cbind(VAR_results, NNS.ARMA_results, `Auto ARMA RMSEs` = auto.ARMA_RMSEs, 
    `Auto ARMA Cor` = auto.ARMA_Correlations), digits = 4) %>% 
    kable_styling(full_width = T) %>% column_spec(1, width = "8cm")

	VAR Estimates RMSE to Actual Values	VAR Estimates Correlation to Actual Values	NNS ARMA RMSEs	NNS ARMA Corr	Auto ARMA RMSEs	Auto ARMA Cor
Diff.Real.GDP	0.5062	-0.3916	0.4194	0.2657	0.464	-0.1329
Interest.Rate	1.3451	0.7741	0.9623	-0.0315	2.631	-0.4812
Inflation	0.3481	-0.6294	0.3255	0.6364	0.227	0.3287

for (i in 1:3) {
    plot(1:h, test_DVs[, i], type = "l", lwd = 3, ylim = c(min(new_values[, 
        i]), max(new_values[, i])), xlab = "Forecast Period", 
        ylab = colnames(new_values)[i])
    lines(1:h, nns_IVs[, i], col = "red", lwd = 2)
    lines(1:h, VAR_test$fcst[[i]][, 1], col = "blue", lwd = 2)
    lines(1:h, arma.fit[[i]], col = "brown", lwd = 2)
    legend("topleft", col = c("black", "red", "blue", "brown"), 
        legend = c("Actual ", "NNS", "VAR", "ARMA"), lty = 1, 
        lwd = c(3, 2, 2, 2))
}

# Significant Difference in Distributions
for (i in 1:3) {
    print(paste0("NNS KS Test: ", suppressWarnings(ks.test(nns_IVs[, 
        i], test_DVs[, i])$p.value)))
    print(paste0("ARMA KS Test: ", suppressWarnings(ks.test(arma.fit[[i]], 
        test_DVs[, i])$p.value)))
}

## [1] "NNS KS Test: 0.868981671175775"
## [1] "ARMA KS Test: 0.0995467717099161"
## [1] "NNS KS Test: 0.0995618483147803"
## [1] "ARMA KS Test: 1.22884247066857e-05"
## [1] "NNS KS Test: 0.00785901405096467"
## [1] "ARMA KS Test: 0.0995467717099161"

7 Use Forecasted IVs in `NNS` and Random Forest

We will now utilize the forecasted variables in the multi-variate NNS.reg() function. Our first step is to pass these variables through the NNS.boost() function in order to reduce any unnecessary features. The objective function obj.fn and objective objective parameters in NNS.boost() need to be re-defined from their default values since this is not a classification problem.

From there, we use the reduced variable set (if any reduction) in the NNS.stack() function which cross- validates the n.best and threshold parameters from the NNS.reg() function.

y_train_nns = head(lagged_new_values, dim(lagged_new_values)[1] - 
    h)
y_test_nns = tail(lagged_new_values, h)

# Select tau = 0 as test set DVs
DVs = which(grepl("tau.0", colnames(y_train_nns)))

nns_est = list()

NNS_est_RMSEs = numeric()
NNS_est_Correlations = numeric()

RF_RMSEs = numeric()
RF_Correlations = numeric()


for (i in DVs) {
    index = which(DVs == i)
    test.set = test_DVs[, index]
    
    # NNS.boost() is an ensemble method comparable to xgboost,
    # and aids in dimension reduction
    nns_boost_est = NNS.boost(IVs.train = y_train_nns[, -i], 
        DV.train = y_train_nns[, i], IVs.test = y_test_nns[, 
            -i], obj.fn = expression(sum((predicted - actual)^2)), 
        objective = "min", ts.test = 2 * h, learner.trials = 100, 
        epochs = 100, ncores = 1, type = NULL, feature.importance = FALSE, 
        folds = 1)
    
    # NNS.stack() cross-validates the parameters of the
    # multivariate NNS.reg() and dimension reduction NNS.reg()
    relevant_vars = colnames(y_train_nns) %in% names(nns_boost_est$feature.weights)
    
    nns_reg_est = NNS.stack(IVs.train = y_train_nns[, relevant_vars], 
        DV.train = y_train_nns[, i], IVs.test = y_test_nns[, 
            relevant_vars], folds = 1, ts.test = 2 * h, order = "max", 
        obj.fn = expression(sum((predicted - actual)^2)), objective = "min")$stack
    
    # Ensemble with univariate estimates
    nns_est[[index]] = (nns_IVs[, index] + nns_reg_est)/2
    
    # Random Forest
    rf = randomForest(y_train_nns[, -i], y_train_nns[, i], ntree = 100)
    rf.pred = predict(rf, newdata = y_test_nns[, -i])
    
    # Print all results
    NNS_est_RMSEs[index] = sqrt(mean((nns_est[[index]] - test.set)^2))
    RF_RMSEs[index] = sqrt(mean((rf.pred - test.set)^2))
    
    NNS_est_Correlations[index] = cor(nns_est[[index]], test.set, 
        method = "spearman")
    RF_Correlations[index] = cor(rf.pred, test.set, method = "spearman")
    
    
    # Plot all results
    plot(1:12, test_DVs[, index], type = "l", lwd = 3, ylim = c(min(new_values[, 
        index]), max(new_values[, index])), xlab = "Forecast Period", 
        ylab = colnames(new_values)[i])
    lines(1:12, nns_est[[index]], col = "red", lwd = 2)
    lines(1:12, rf.pred, col = "green", lwd = 2)
    lines(1:12, VAR_test$fcst[[index]][, 1], col = "blue", lwd = 2)
    legend("topleft", col = c("black", "red", "blue", "green"), 
        legend = c("Actual ", "NNS", "VAR", "RF"), lty = 1, lwd = c(3, 
            2, 2, 2))
}

ALL_RESULTS = cbind(VAR_results, NNS.ARMA_results, `NNS RMSEs` = NNS_est_RMSEs, 
    `NNS Corr` = NNS_est_Correlations, `RF RMSEs` = RF_RMSEs, 
    `RF Corr` = RF_Correlations)

knitr::kable(ALL_RESULTS, digits = 4) %>% kable_styling(full_width = T) %>% 
    column_spec(1, width = "8cm")

	VAR Estimates RMSE to Actual Values	VAR Estimates Correlation to Actual Values	NNS ARMA RMSEs	NNS ARMA Corr	NNS RMSEs	NNS Corr	RF RMSEs	RF Corr
Diff.Real.GDP	0.5062	-0.3916	0.4194	0.2657	0.4335	0.3007	0.4949	-0.3497
Interest.Rate	1.3451	0.7741	0.9623	-0.0315	0.6032	0.5814	0.5211	0.7110
Inflation	0.3481	-0.6294	0.3255	0.6364	0.3031	0.5804	0.2631	0.4056

8 `NNS.VAR()`

The NNS.VAR() function now accomplishes all of the prior steps in a single line of code. NNS.VAR does take the additional step of weighting each estimate by its objective function result.

nns_var_estimates = NNS.VAR(y_train_VAR, h = 12, tau = 4, ncores = 1)

nns_var_estimates

## $univariate
##       Diff Real GDP Interest Rate Inflation
##  [1,]    0.71182079      4.536222 0.6891101
##  [2,]    0.97198238      4.655889 0.5132415
##  [3,]    0.49344043      4.603333 0.5697954
##  [4,]    0.89140137      4.737889 0.4357171
##  [5,]    1.00169992      4.871167 0.3667949
##  [6,]    0.93692811      4.665000 0.4023420
##  [7,]    0.04796763      4.750889 0.4080647
##  [8,]    0.56628191      4.603222 0.4370086
##  [9,]    0.41366022      4.474333 0.4392642
## [10,]    0.56131736      4.821333 0.3732100
## [11,]    0.67121142      4.492222 0.3300415
## [12,]    0.77318203      4.568333 0.2547801
## 
## $multivariate
##       Diff Real GDP Interest Rate Inflation
##  [1,]     0.9279022      2.358849 0.4422726
##  [2,]     0.9141970      2.440993 0.5782729
##  [3,]     0.9844185      4.717534 0.7800222
##  [4,]     0.8247261      4.666915 0.5522493
##  [5,]     0.8331298      4.675716 0.4952923
##  [6,]     0.9397740      5.092680 0.4653352
##  [7,]     1.1505812      5.161801 0.5074708
##  [8,]     0.9506566      4.965832 0.3940349
##  [9,]     1.0781853      4.825252 0.4178213
## [10,]     1.1659552      5.145660 0.5238891
## [11,]     1.2304537      4.768645 0.3865078
## [12,]     0.9972291      4.856568 0.3262525
## 
## $ensemble
##       Diff Real GDP Interest Rate Inflation
##  [1,]     0.7613002      2.759611 0.5048272
##  [2,]     0.9587504      2.848662 0.5617924
##  [3,]     0.6058670      4.696514 0.7267456
##  [4,]     0.8761337      4.679979 0.5227172
##  [5,]     0.9630999      4.711690 0.4627279
##  [6,]     0.9375798      5.013962 0.4493712
##  [7,]     0.3004495      5.086169 0.4822789
##  [8,]     0.6542979      4.899091 0.4049255
##  [9,]     0.5658265      4.760663 0.4232554
## [10,]     0.6997703      5.085965 0.4857034
## [11,]     0.7992695      4.717768 0.3721978
## [12,]     0.8244854      4.803516 0.3081397

NNS_VAR_RMSEs = numeric()
NNS_VAR_corr = numeric()

for (i in 1:3) {
    NNS_VAR_RMSEs[i] = sqrt(mean((nns_var_estimates$ensemble[, 
        i] - test_DVs[, i])^2))
    
    NNS_VAR_corr[i] = cor(nns_var_estimates$ensemble[, i], test_DVs[, 
        i], method = "spearman")
}

NNS_VAR_results = cbind(NNS_VAR_RMSEs, NNS_VAR_corr)
rownames(NNS_VAR_results) = colnames(y_train_VAR)
knitr::kable(NNS_VAR_results, digits = 4) %>% kable_styling(full_width = T) %>% 
    column_spec(1, width = "8cm")

	NNS_VAR_RMSEs	NNS_VAR_corr
Diff Real GDP	0.4082	0.2657
Interest Rate	0.4855	0.8231
Inflation	0.2963	0.5175

Data available for download at the following: https://sites.google.com/site/lkilian2019/↩︎

Multivariate Time Series Forecasting: Nonparametric Vector Autoregression

Fred Viole

1 Introduction

2 `NNS.VAR` Example with 3 Variables: Fed Funds, Real GNP, and Inflation

3 Forecast 12 Quarters Out-of-Sample Using All Methods

4 Traditional VAR

5 `NNS` and Random Forest

5.1 Function to Generate Lagged Variables

6 Predictions

6.1 1 Year Lag \((\tau = 4)\), 3 Year Forecast \((h = 12)\)

6.2 Create Forecasted IVs Using `NNS.ARMA()`

7 Use Forecasted IVs in `NNS` and Random Forest

8 `NNS.VAR()`

Multivariate Time Series Forecasting: Nonparametric Vector Autoregression

Fred Viole

1 Introduction

2 NNS.VAR Example with 3 Variables: Fed Funds, Real GNP, and Inflation

3 Forecast 12 Quarters Out-of-Sample Using All Methods

4 Traditional VAR

5 NNS and Random Forest

5.1 Function to Generate Lagged Variables

6 Predictions

6.1 1 Year Lag \((\tau = 4)\), 3 Year Forecast \((h = 12)\)

6.2 Create Forecasted IVs Using NNS.ARMA()

7 Use Forecasted IVs in NNS and Random Forest

8 NNS.VAR()

2 `NNS.VAR` Example with 3 Variables: Fed Funds, Real GNP, and Inflation

5 `NNS` and Random Forest

6.2 Create Forecasted IVs Using `NNS.ARMA()`

7 Use Forecasted IVs in `NNS` and Random Forest

8 `NNS.VAR()`