library("midasr") # This script replicates results obtained from the MIDAS # Matlab Toolbox quite faithfully. It is based on code # taken from the midasr documentation, but the sample range # and lag specification are adjusted so as to match Matlab. # Ghysels' data (taken from mydata.xlsx) load("qgdp.Rdata") load("payems.Rdata") USqgdp <- qgdp USpayems <- payems y <- window(USqgdp, start = c(1947, 1), end = c(2011, 2)) x <- window(USpayems, start = c(1947, 1), end = c(2011, 7)) yg <- diff(log(y)) * 100 xg <- diff(log(x)) * 100 ny <- ts(yg, start=c(1947, 2), frequency=4) nx <- ts(xg, start=c(1947, 2), frequency=12) yy <- window(ny, start=c(1984,1), end=c(2009,1)) xx <- window(nx, start=c(1984,1), end=c(2009,3)) cat("\n=== beta0 ===\n") beta0 <- midas_r(yy ~ mls(yy, 1, 1) + mls(xx, 5:13, 3, nbeta), start = list(xx = c(1.7, 1, 5))) summary(beta0) uhat <- beta0$residuals sprintf("T = %d, SSR = %.8f", beta0$nobs, sum(uhat^2)) cat("\n=== betaN ===\n") betan <- midas_r(yy ~ mls(yy, 1, 1) + mls(xx, 5:13, 3, nbetaMT), start = list(xx = c(2, 1, 5, 0))) summary(betan) uhat <- betan$residuals sprintf("T = %d, SSR = %.8f", betan$nobs, sum(uhat^2)) cat("\n=== U-MIDAS ===\n") um <- midas_r(yy ~ mls(yy, 1, 1) + mls(xx, 5:13, 3), start = NULL) summary(um) uhat <- um$residuals sprintf("T = %d, SSR = %.8f", um$nobs, sum(uhat^2)) # RMSEs via average_forecast() fulldata <- list(xx = window(nx, start=c(1984,1), end=c(2011,6)), yy = window(ny, start=c(1984,1), end=c(2011,2))) insample <- 1:length(yy) outsample <- (1:length(fulldata$yy))[-insample] avgf <- average_forecast(list(beta0, betan, um), fulldata, insample, outsample) sqrt(avgf$accuracy$individual$MSE.out.of.sample)