library(quantmod)
## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## Loading required package: TTR
## Version 0.4-0 included new data defaults. See ?getSymbols.
gdp <- getSymbols("NAEXKP01CAQ189S", src = "FRED", auto.assign = FALSE)
## As of 0.4-0, 'getSymbols' uses env=parent.frame() and
## auto.assign=TRUE by default.
##
## This behavior will be phased out in 0.5-0 when the call will
## default to use auto.assign=FALSE. getOption("getSymbols.env") and
## getOptions("getSymbols.auto.assign") are now checked for alternate defaults
##
## This message is shown once per session and may be disabled by setting
## options("getSymbols.warning4.0"=FALSE). See ?getSymbol for more details
chartSeries(gdp ,theme="white", name = "Total GDP for Canada")
# Y GDP available till 3rd quarter in 2014
y <- window(ts(gdp, start = c(1981, 1), frequency = 4), end = c(2014,3))
# log differenced
yg <- log(y/lag(y, 1)) * 100
#equalize the sample sizes, which are different in the original
#data. simply add additional NA values at the beginning and the end of the data.
ny <- ts(c(NA,yg,NA), start = start(y), frequency = 4)
# might use na.locf.xts {xts} xts method to replace ‘NA’ with most recent non-‘NA’
# training data from 1981 to 2011
yy <- window(ny, start = c(1981, 1), end = c(2011, 4))
#length(yy) #124
labor <- getSymbols("LFEMTTTTCAM647N", src = "FRED", auto.assign = FALSE)
chartSeries(labor ,theme="white", name = "Employed Population for Canada")
#
x <- window(ts(labor, start = c(1955, 1), frequency = 12), end = c(2014,12))
# differenced
xg <- log(x/lag(x, 1)) * 100
# equalize the sample sizes with
nx <- ts(c(NA,xg), start = start(x), frequency = 12)
#length(xx) # 372/12
# training data
xx <- window(nx, start = c(1981, 1), end = c(2011, 12))
tsx <- getSymbols("^GSPTSE",src="yahoo", from = '1977-01-01', auto.assign = FALSE)
chartSeries(tsx ,theme="white", name = "S&P/TSX Composite index")
difftsx <- log(tsx/lag(tsx, 1)) * 100
library(zoo)
tsx.zoo <- zoo(difftsx$GSPTSE.Adjusted)
# define a function to pad NAs at the beginning of the each month.
padd_nas <- function(x, desired_length) {
n <- length(x)
if(n < desired_length) {
c(rep(NA,desired_length-n),x)
} else {
tail(x,desired_length)
}
}
# take 22 days entries for each month
#tsxdaily <- aggregate(tsx.zoo, as.yearmon, padd_nas, 22)
# zoo method as.yearqtr working
tsxdaily <- aggregate(tsx.zoo, as.yearqtr, padd_nas, 66)
# The zoo package has a yearqtr class that can represent quarterly
#data. It is also possible to create zoo series of quarterly data. Such
#series can be converted to ts class and that class can some
#builtin support for quarterly data as well.
# tsxdaily is zoo series index by yearmon or yearqtr, so it is like a matrix. we have to vectorize it with t() transpose function since vector function take element by column.
#zz <- as.vector(t(window(tsxdaily, start = as.yearmon("1981-1"),
# end =as.yearmon("2011-12"))))
## for a quarterly object, this is more proper.
# training data
zz <- as.vector(t(window(tsxdaily, start = as.yearqtr("1981 Q1"),
end =as.yearqtr("2011 Q4"))))
The unrestricted MIDAS regression is similar to OLS without weight scheme.
\[{y_t= y_{t-1} + \sum_{h=0}^d\theta_{h}x_{tm-h}+\sum_{j=0}^p\phi_{j}z_{tn-j}+u_t}\]
library("midasr")
## Loading required package: sandwich
## Loading required package: optimx
# only monthly labor
umt <- midas_r(yy ~ mls(yy, 1, 1) + fmls(xx, 5, 3) , start = NULL)
coef(umt)
## (Intercept) yy xx1 xx2 xx3 xx4
## -0.20577481 0.37353457 0.40204349 0.20159802 -0.06639832 -0.03751980
## xx5 xx6
## 0.27419057 0.59385162
Please see appendix 1 for the lag selection of monthly labor variable.
normalized exponential Almon lag restricts the coefficients theta_h in the following way:
\[{\theta_{h}=\delta\frac{\exp(\lambda_1(h+1)+\dots+\lambda_r(h+1)^r)}{\sum_{s=0}^d\exp(\lambda_1(s+1)+\dots+\lambda_r(h+1)^r)}}\]
The parameter d should be the first element in vector p. The degree of the polynomial is then decided by the number of the remaining parameters.
trend <- 1:length(yy)
nealmont <- midas_r(midas_r(yy~mls(yy, 1, 1)+fmls(xx,11,3,nealmon),
start=list(xx=rep(0,3))))
coef(nealmont)
## (Intercept) yy xx1 xx2 xx3
## -0.209800428 0.345094033 1.669145445 -0.037013549 0.002411721
beta0t <- midas_r(yy ~ mls(yy, 1, 1) + mls(xx, 0:5, 3, nbeta),
start = list(xx = c(1.7, 1, 5)))
#beta0 <- midas_r(beta0,Ofunction="nls")
coef(beta0t)
## (Intercept) yy xx1 xx2 xx3
## -0.25333044 0.56315435 -0.05458931 16.53215824 24.56574949
# Test if the variable is adequate
#hAh.test(beta0)
betant <- midas_r(yy ~ mls(yy, 1, 1) + mls(xx, 0:11, 3, nbetaMT),
start = list(xx = c(2, 1, 5, 0)))
coef(betant)
## (Intercept) yy xx1 xx2 xx3 xx4
## -0.2087275 0.3537060 1.6270627 24.3028280 4.9067814 52.9835027
#hAh.test(betant)
library(zoo)
fulldata <- list(xx=window(nx,start=c(1981,1),end=c(2013,12)),
yy=window(ny,start=c(1981,1),end=c(2013,4)),
zz = as.vector(t(window(tsxdaily,
start = as.yearqtr("1981-1"),
end =as.yearqtr("2013-4")))))
# training sample
insample <- 1:length(yy)
# testing sample
outsample <- (1:length(fulldata$yy))[-insample]
avgf<-average_forecast( list(umt, beta0t,betant, nealmont),
data=fulldata,insample=insample,outsample=outsample)
sqrt(avgf$accuracy$individual$MSE.out.of.sample)
## [1] 0.5944451 0.2622549 0.2558328 0.2575742
Beta with non-zero weight specifications weight specifications is the best.
# yy, xx, zz, only 22 days daily data can be used due to the limit of degree of freedom
umd <- midas_r(yy ~ mls(yy, 1, 1) + fmls(xx, 11, 3) +
fmls(zz, 21, 66) , start = NULL)
coef(umd)
## (Intercept) yy xx1 xx2 xx3
## -0.207883714 0.304816354 0.510311600 0.361960658 -0.004820447
## xx4 xx5 xx6 xx7 xx8
## 0.090774437 0.314276919 0.585058482 0.012122081 -0.099216438
## xx9 xx10 xx11 xx12 zz1
## -0.033239355 0.256912261 -0.022005467 -0.382471926 0.094292498
## zz2 zz3 zz4 zz5 zz6
## 0.061209686 0.027910026 0.151505832 0.084230389 -0.054593787
## zz7 zz8 zz9 zz10 zz11
## 0.118138031 -0.004954854 0.087646158 0.017288756 -0.098643853
## zz12 zz13 zz14 zz15 zz16
## 0.014917448 -0.036306535 0.000739204 0.009395796 -0.071708096
## zz17 zz18 zz19 zz20 zz21
## 0.015488527 0.005902869 -0.059555486 -0.056478668 0.010923438
## zz22
## -0.150557405
Not enough degree of freedom to run more than 22 lag for daily highfrequency data. (one month has 22 trading days )
nealmond <- midas_r(midas_r(yy~mls(yy, 1, 1)+fmls(xx,11,3,nealmon) +
fmls(zz, 21, 66), start=list(xx=rep(0,3))))
coef(nealmont)
## (Intercept) yy xx1 xx2 xx3
## -0.209800428 0.345094033 1.669145445 -0.037013549 0.002411721
beta0d <- midas_r(yy ~ mls(yy, 1, 1) + mls(xx, 0:5, 3, nbeta) +
fmls(zz, 8, 66, nbeta),
start = list(xx = c(1.7, 1, 5),
zz = c(1.7, 1, 5)))
coef(beta0d)
## (Intercept) yy xx1 xx2 xx3 zz1
## -0.26379636 0.55724133 0.04359155 6.65608463 -1.68459265 0.13497652
## zz2 zz3
## 11.50151312 37.72330642
# Test if the variable is adequate
#hAh.test(beta0)
betand <- midas_r(yy ~ mls(yy, 1, 1) + fmls(xx, 11, 3, nbetaMT) +
fmls(zz, 34, 66, nbetaMT),
start = list(xx = c(2, 1, 5, 0), zz = c(2, 1, 5, 0)))
coef(betand)
## (Intercept) yy xx1 xx2 xx3 xx4
## -0.2062569 0.3334300 1.6341859 21.0440597 -18.5421236 51.1414737
## zz1 zz2 zz3 zz4
## -0.7449314 -2.6980656 4.9957138 26.4169209
#hAh.test(betant)
avgfd<-average_forecast( list(umd, nealmond, beta0d,betand),
data=fulldata,insample=insample,outsample=outsample)
sqrt(avgfd$accuracy$individual$MSE.out.of.sample)
## [1] 0.8578126 0.3720912 0.2934550 0.2619594
Beta with non-zero weight specifications is the best.
# Table gives AIC BIC to pick up a suitable lag structure
trend <- 1:1:length(yy)
mlr <- hf_lags_table(yy ~ mls(yy, 1, 1) + mls(xx, 0:5, 3, nbeta),
start = list(xx = c(1.7, 1, 5)),
from=c(xx=0),to=list(xx=c(3,12)))
##
## Model selection progress:
##
|
| | 0%
|
|====== | 10%
|
|============= | 20%
|
|==================== | 30%
|
|========================== | 40%
|
|================================ | 50%
|
|======================================= | 60%
|
|============================================== | 70%
|
|==================================================== | 80%
|
|========================================================== | 90%
|
|=================================================================| 100%
mlr
## model AIC.restricted
## 1 yy ~ mls(yy, 1, 1) + mls(xx, 0:3, 3, nbeta) 237.9921
## 2 yy ~ mls(yy, 1, 1) + mls(xx, 0:4, 3, nbeta) 237.7234
## 3 yy ~ mls(yy, 1, 1) + mls(xx, 0:5, 3, nbeta) 237.1546
## 4 yy ~ mls(yy, 1, 1) + mls(xx, 0:6, 3, nbeta) 237.1724
## 5 yy ~ mls(yy, 1, 1) + mls(xx, 0:7, 3, nbeta) 237.1403
## 6 yy ~ mls(yy, 1, 1) + mls(xx, 0:8, 3, nbeta) 237.8306
## 7 yy ~ mls(yy, 1, 1) + mls(xx, 0:9, 3, nbeta) 237.6873
## 8 yy ~ mls(yy, 1, 1) + mls(xx, 0:10, 3, nbeta) 237.7061
## 9 yy ~ mls(yy, 1, 1) + mls(xx, 0:11, 3, nbeta) 237.1425
## 10 yy ~ mls(yy, 1, 1) + mls(xx, 0:12, 3, nbeta) 237.5524
## BIC.restricted AIC.unrestricted BIC.unrestricted hAh.test.p.value First
## 1 254.7171 238.6546 258.1671 2.627309e-01 FALSE
## 2 254.4483 234.7382 257.0381 4.403120e-02 FALSE
## 3 253.8796 220.5251 245.6125 4.313015e-05 FALSE
## 4 253.8974 222.1830 250.0579 1.034661e-04 FALSE
## 5 253.8652 222.3204 252.9829 1.232867e-04 FALSE
## 6 254.5555 223.6811 257.1310 1.768001e-04 FALSE
## 7 254.4123 217.5726 253.8100 8.974484e-06 FALSE
## 8 254.4310 217.1652 256.1901 7.254784e-06 FALSE
## 9 253.8675 211.2996 253.1120 6.613745e-07 FALSE
## 10 254.2773 206.8370 251.4369 6.915395e-08 FALSE
## Second Convergence
## 1 FALSE 0
## 2 FALSE 0
## 3 TRUE 0
## 4 FALSE 0
## 5 TRUE 0
## 6 TRUE 0
## 7 FALSE 0
## 8 TRUE 0
## 9 TRUE 1
## 10 TRUE 0
# 5 lag is the best
mlr <- hf_lags_table(yy ~ mls(yy, 1, 1) + mls(xx, 0:11, 3, nbetaMT),
start = list(xx = c(2, 1, 5, 0)),
from=c(xx=0),to=list(xx=c(3,12)))
##
## Model selection progress:
##
|
| | 0%
|
|======= | 11%
|
|============== | 22%
|
|====================== | 33%
|
|============================= | 44%
|
|==================================== | 56%
|
|=========================================== | 67%
|
|=================================================== | 78%
|
|========================================================== | 89%
|
|=================================================================| 100%
mlr
## model AIC.restricted
## 1 yy ~ mls(yy, 1, 1) + mls(xx, 0:4, 3, nbetaMT) 240.0955
## 2 yy ~ mls(yy, 1, 1) + mls(xx, 0:5, 3, nbetaMT) 239.4388
## 3 yy ~ mls(yy, 1, 1) + mls(xx, 0:6, 3, nbetaMT) 240.1737
## 4 yy ~ mls(yy, 1, 1) + mls(xx, 0:7, 3, nbetaMT) 240.3205
## 5 yy ~ mls(yy, 1, 1) + mls(xx, 0:8, 3, nbetaMT) 240.3215
## 6 yy ~ mls(yy, 1, 1) + mls(xx, 0:9, 3, nbetaMT) 238.8796
## 7 yy ~ mls(yy, 1, 1) + mls(xx, 0:10, 3, nbetaMT) 236.0066
## 8 yy ~ mls(yy, 1, 1) + mls(xx, 0:11, 3, nbetaMT) 231.5997
## 9 yy ~ mls(yy, 1, 1) + mls(xx, 0:12, 3, nbetaMT) 237.6039
## BIC.restricted AIC.unrestricted BIC.unrestricted hAh.test.p.value First
## 1 259.6079 234.7382 257.0381 7.517600e-03 FALSE
## 2 258.9513 220.5251 245.6125 7.639384e-06 FALSE
## 3 259.6862 222.1830 250.0579 1.903407e-05 FALSE
## 4 259.8330 222.3204 252.9829 2.380760e-05 FALSE
## 5 259.8339 223.6811 257.1310 5.467211e-05 FALSE
## 6 258.3920 217.5726 253.8100 5.267758e-06 FALSE
## 7 255.5190 217.1652 256.1901 1.974658e-05 FALSE
## 8 251.1121 211.2996 253.1120 9.557245e-06 FALSE
## 9 257.1163 206.8370 251.4369 3.412782e-08 FALSE
## Second Convergence
## 1 FALSE 0
## 2 FALSE 0
## 3 FALSE 0
## 4 FALSE 0
## 5 FALSE 0
## 6 FALSE 0
## 7 FALSE 0
## 8 FALSE 0
## 9 FALSE 0
# 11 lag is the best
mlr <- hf_lags_table(yy~trend+fmls(xx,11,3,nealmon),
start=list(xx=rep(0,3)),
from=c(xx=0),to=list(xx=c(2,12)))
##
## Model selection progress:
##
|
| | 0%
|
|====== | 10%
|
|============= | 20%
|
|==================== | 30%
|
|========================== | 40%
|
|================================ | 50%
|
|======================================= | 60%
|
|============================================== | 70%
|
|==================================================== | 80%
|
|========================================================== | 90%
|
|=================================================================| 100%
mlr
## model AIC.restricted BIC.restricted
## 1 yy ~ trend + mls(xx, 0:3, 3, nealmon) 281.4458 298.1708
## 2 yy ~ trend + mls(xx, 0:4, 3, nealmon) 278.1852 294.9101
## 3 yy ~ trend + mls(xx, 0:5, 3, nealmon) 239.7061 256.4311
## 4 yy ~ trend + mls(xx, 0:6, 3, nealmon) 278.8665 295.5915
## 5 yy ~ trend + mls(xx, 0:7, 3, nealmon) 277.5901 294.3150
## 6 yy ~ trend + mls(xx, 0:8, 3, nealmon) 275.7834 292.5084
## 7 yy ~ trend + mls(xx, 0:9, 3, nealmon) 253.6401 270.3651
## 8 yy ~ trend + mls(xx, 0:10, 3, nealmon) 243.5855 260.3105
## 9 yy ~ trend + mls(xx, 0:11, 3, nealmon) 238.5657 255.2907
## 10 yy ~ trend + mls(xx, 0:12, 3, nealmon) 246.9708 263.6958
## AIC.unrestricted BIC.unrestricted hAh.test.p.value First Second
## 1 283.3236 302.8360 7.339045e-01 FALSE TRUE
## 2 274.5653 296.8652 2.462407e-02 FALSE TRUE
## 3 243.2915 268.3789 5.170606e-01 FALSE TRUE
## 4 245.1850 273.0599 2.348721e-09 FALSE TRUE
## 5 236.3078 266.9702 2.308487e-11 FALSE TRUE
## 6 234.0998 267.5497 2.234357e-11 FALSE TRUE
## 7 222.5792 258.8166 2.056697e-08 FALSE TRUE
## 8 223.4070 262.4319 8.745710e-06 FALSE TRUE
## 9 216.3536 258.1660 3.026370e-06 FALSE TRUE
## 10 214.4491 259.0490 1.004312e-08 FALSE TRUE
## Convergence
## 1 0
## 2 0
## 3 0
## 4 0
## 5 0
## 6 0
## 7 0
## 8 0
## 9 0
## 10 0
# 11 las is the best
#-------------------
# For zz
# Table gives AIC BIC to pick up a suitable lag structure
trend <- 1:1:length(yy)
mlr <- hf_lags_table(yy ~ mls(yy, 1, 1) + mls(zz, 0:5, 66, nbeta),
start = list(zz = c(1.7, 1, 5)),
from=c(zz=0),to=list(zz=c(5,56)))
##
## Model selection progress:
##
|
| | 0%
|
|= | 2%
|
|== | 4%
|
|==== | 6%
|
|===== | 8%
|
|====== | 10%
|
|======== | 12%
|
|========= | 13%
|
|========== | 15%
|
|=========== | 17%
|
|============ | 19%
|
|============== | 21%
|
|=============== | 23%
|
|================ | 25%
|
|================== | 27%
|
|=================== | 29%
|
|==================== | 31%
|
|===================== | 33%
|
|====================== | 35%
|
|======================== | 37%
|
|========================= | 38%
|
|========================== | 40%
|
|============================ | 42%
|
|============================= | 44%
|
|============================== | 46%
|
|=============================== | 48%
|
|================================ | 50%
|
|================================== | 52%
|
|=================================== | 54%
|
|==================================== | 56%
|
|====================================== | 58%
|
|======================================= | 60%
|
|======================================== | 62%
|
|========================================= | 63%
|
|========================================== | 65%
|
|============================================ | 67%
|
|============================================= | 69%
|
|============================================== | 71%
|
|================================================ | 73%
|
|================================================= | 75%
|
|================================================== | 77%
|
|=================================================== | 79%
|
|==================================================== | 81%
|
|====================================================== | 83%
|
|======================================================= | 85%
|
|======================================================== | 87%
|
|========================================================== | 88%
|
|=========================================================== | 90%
|
|============================================================ | 92%
|
|============================================================= | 94%
|
|============================================================== | 96%
|
|================================================================ | 98%
|
|=================================================================| 100%
# 8 lag is the best
mlr
## model AIC.restricted
## 1 yy ~ mls(yy, 1, 1) + mls(zz, 0:5, 66, nbeta) 241.3531
## 2 yy ~ mls(yy, 1, 1) + mls(zz, 0:6, 66, nbeta) 244.4327
## 3 yy ~ mls(yy, 1, 1) + mls(zz, 0:7, 66, nbeta) 244.4378
## 4 yy ~ mls(yy, 1, 1) + mls(zz, 0:8, 66, nbeta) 240.9166
## 5 yy ~ mls(yy, 1, 1) + mls(zz, 0:9, 66, nbeta) 241.2049
## 6 yy ~ mls(yy, 1, 1) + mls(zz, 0:10, 66, nbeta) 241.2420
## 7 yy ~ mls(yy, 1, 1) + mls(zz, 0:11, 66, nbeta) 241.2571
## 8 yy ~ mls(yy, 1, 1) + mls(zz, 0:12, 66, nbeta) 241.2691
## 9 yy ~ mls(yy, 1, 1) + mls(zz, 0:13, 66, nbeta) 241.2808
## 10 yy ~ mls(yy, 1, 1) + mls(zz, 0:14, 66, nbeta) 241.2911
## 11 yy ~ mls(yy, 1, 1) + mls(zz, 0:15, 66, nbeta) 241.3001
## 12 yy ~ mls(yy, 1, 1) + mls(zz, 0:16, 66, nbeta) 241.3077
## 13 yy ~ mls(yy, 1, 1) + mls(zz, 0:17, 66, nbeta) 241.3143
## 14 yy ~ mls(yy, 1, 1) + mls(zz, 0:18, 66, nbeta) 241.3201
## 15 yy ~ mls(yy, 1, 1) + mls(zz, 0:19, 66, nbeta) 241.3251
## 16 yy ~ mls(yy, 1, 1) + mls(zz, 0:20, 66, nbeta) 241.3295
## 17 yy ~ mls(yy, 1, 1) + mls(zz, 0:21, 66, nbeta) 241.3335
## 18 yy ~ mls(yy, 1, 1) + mls(zz, 0:22, 66, nbeta) 241.3371
## 19 yy ~ mls(yy, 1, 1) + mls(zz, 0:23, 66, nbeta) 241.3404
## 20 yy ~ mls(yy, 1, 1) + mls(zz, 0:24, 66, nbeta) 241.3433
## 21 yy ~ mls(yy, 1, 1) + mls(zz, 0:25, 66, nbeta) 241.3460
## 22 yy ~ mls(yy, 1, 1) + mls(zz, 0:26, 66, nbeta) 241.3485
## 23 yy ~ mls(yy, 1, 1) + mls(zz, 0:27, 66, nbeta) 241.3508
## 24 yy ~ mls(yy, 1, 1) + mls(zz, 0:28, 66, nbeta) 241.3529
## 25 yy ~ mls(yy, 1, 1) + mls(zz, 0:29, 66, nbeta) 241.3549
## 26 yy ~ mls(yy, 1, 1) + mls(zz, 0:30, 66, nbeta) 241.3568
## 27 yy ~ mls(yy, 1, 1) + mls(zz, 0:31, 66, nbeta) 241.3585
## 28 yy ~ mls(yy, 1, 1) + mls(zz, 0:32, 66, nbeta) 241.3601
## 29 yy ~ mls(yy, 1, 1) + mls(zz, 0:33, 66, nbeta) 241.3616
## 30 yy ~ mls(yy, 1, 1) + mls(zz, 0:34, 66, nbeta) 241.3632
## 31 yy ~ mls(yy, 1, 1) + mls(zz, 0:35, 66, nbeta) 241.3648
## 32 yy ~ mls(yy, 1, 1) + mls(zz, 0:36, 66, nbeta) 241.3666
## 33 yy ~ mls(yy, 1, 1) + mls(zz, 0:37, 66, nbeta) 241.3672
## 34 yy ~ mls(yy, 1, 1) + mls(zz, 0:38, 66, nbeta) 241.3693
## 35 yy ~ mls(yy, 1, 1) + mls(zz, 0:39, 66, nbeta) 241.3692
## 36 yy ~ mls(yy, 1, 1) + mls(zz, 0:40, 66, nbeta) 241.3701
## 37 yy ~ mls(yy, 1, 1) + mls(zz, 0:41, 66, nbeta) 241.3711
## 38 yy ~ mls(yy, 1, 1) + mls(zz, 0:42, 66, nbeta) 241.3749
## 39 yy ~ mls(yy, 1, 1) + mls(zz, 0:43, 66, nbeta) 241.3729
## 40 yy ~ mls(yy, 1, 1) + mls(zz, 0:44, 66, nbeta) 241.3737
## 41 yy ~ mls(yy, 1, 1) + mls(zz, 0:45, 66, nbeta) 241.3745
## 42 yy ~ mls(yy, 1, 1) + mls(zz, 0:46, 66, nbeta) 241.3753
## 43 yy ~ mls(yy, 1, 1) + mls(zz, 0:47, 66, nbeta) 241.3768
## 44 yy ~ mls(yy, 1, 1) + mls(zz, 0:48, 66, nbeta) 241.3768
## 45 yy ~ mls(yy, 1, 1) + mls(zz, 0:49, 66, nbeta) 241.3774
## 46 yy ~ mls(yy, 1, 1) + mls(zz, 0:50, 66, nbeta) 241.3781
## 47 yy ~ mls(yy, 1, 1) + mls(zz, 0:51, 66, nbeta) 241.3787
## 48 yy ~ mls(yy, 1, 1) + mls(zz, 0:52, 66, nbeta) 241.3795
## 49 yy ~ mls(yy, 1, 1) + mls(zz, 0:53, 66, nbeta) 241.3799
## 50 yy ~ mls(yy, 1, 1) + mls(zz, 0:54, 66, nbeta) 241.3807
## 51 yy ~ mls(yy, 1, 1) + mls(zz, 0:55, 66, nbeta) 241.3812
## 52 yy ~ mls(yy, 1, 1) + mls(zz, 0:56, 66, nbeta) 241.3817
## BIC.restricted AIC.unrestricted BIC.unrestricted hAh.test.p.value First
## 1 258.1772 245.0055 270.2417 0.52901277 FALSE
## 2 261.2569 245.4411 273.4813 0.15490081 FALSE
## 3 261.2619 247.2446 278.0889 0.23644450 FALSE
## 4 257.7407 248.0420 281.6903 0.60603115 FALSE
## 5 258.0290 249.0910 285.5433 0.58080608 FALSE
## 6 258.0661 251.0383 290.2946 0.68236522 FALSE
## 7 258.0812 250.7904 292.8507 0.55836671 FALSE
## 8 258.0932 252.6972 297.5616 0.64950081 FALSE
## 9 258.1049 254.4659 302.1342 0.71846144 FALSE
## 10 258.1152 255.6085 306.0809 0.73057616 FALSE
## 11 258.1242 256.4815 309.7579 0.72277958 FALSE
## 12 258.1319 258.4765 314.5569 0.79341362 FALSE
## 13 258.1385 260.2902 319.1747 0.84138946 FALSE
## 14 258.1442 250.8327 312.5212 0.20093694 FALSE
## 15 258.1492 248.8087 313.3012 0.11208691 FALSE
## 16 258.1537 250.0395 317.3360 0.13002527 FALSE
## 17 258.1576 247.2083 317.3088 0.05820553 FALSE
## 18 258.1612 248.7649 321.6695 0.07513804 FALSE
## 19 258.1645 249.0108 324.7194 0.07044003 FALSE
## 20 258.1674 248.8252 327.3378 0.05978736 FALSE
## 21 258.1701 249.2644 330.5811 0.05909748 FALSE
## 22 258.1726 249.7398 333.8604 0.05904890 FALSE
## 23 258.1749 251.5307 338.4554 0.07939698 FALSE
## 24 258.1771 253.5218 343.2505 0.10878345 FALSE
## 25 258.1790 255.0657 347.5984 0.13345191 FALSE
## 26 258.1809 257.0655 352.4022 0.17493281 FALSE
## 27 258.1826 259.0387 357.1794 0.22263887 FALSE
## 28 258.1842 256.6729 357.6177 0.13766019 FALSE
## 29 258.1858 258.6673 362.4160 0.17912656 FALSE
## 30 258.1873 260.1923 366.7451 0.21193593 FALSE
## 31 258.1889 260.2600 369.6168 0.19965112 FALSE
## 32 258.1907 260.0672 372.2281 0.18105906 FALSE
## 33 258.1913 257.6398 372.6046 0.11262184 FALSE
## 34 258.1934 258.0838 375.8527 0.11421502 FALSE
## 35 258.1934 260.0118 380.5847 0.14884173 FALSE
## 36 258.1943 261.7802 385.1571 0.18541583 FALSE
## 37 258.1953 263.7217 389.9027 0.23256311 FALSE
## 38 258.1990 265.7180 394.7030 0.28776249 FALSE
## 39 258.1970 265.7014 397.4904 0.27523476 FALSE
## 40 258.1979 267.6868 402.2798 0.33444304 FALSE
## 41 258.1987 267.8046 405.2017 0.32719727 FALSE
## 42 258.1994 269.7455 409.9466 0.38888869 FALSE
## 43 258.2009 271.6049 414.6100 0.45041458 FALSE
## 44 258.2009 273.3985 419.2076 0.51059667 FALSE
## 45 258.2016 275.1690 423.7821 0.56982128 FALSE
## 46 258.2022 276.6666 428.0837 0.61733911 FALSE
## 47 258.2028 278.4310 432.6521 0.67303553 FALSE
## 48 258.2037 280.3540 437.3792 0.73054413 FALSE
## 49 258.2040 275.5277 435.3569 0.54935081 FALSE
## 50 258.2049 269.8066 432.4398 0.33518351 FALSE
## 51 258.2053 271.6558 437.0931 0.39681100 FALSE
## 52 258.2058 273.6310 441.8722 0.46606143 FALSE
## Second Convergence
## 1 TRUE 0
## 2 TRUE 0
## 3 TRUE 0
## 4 TRUE 0
## 5 TRUE 0
## 6 TRUE 0
## 7 TRUE 0
## 8 TRUE 0
## 9 TRUE 0
## 10 TRUE 0
## 11 TRUE 0
## 12 TRUE 0
## 13 TRUE 0
## 14 TRUE 0
## 15 TRUE 0
## 16 TRUE 0
## 17 TRUE 0
## 18 TRUE 0
## 19 TRUE 0
## 20 TRUE 0
## 21 TRUE 0
## 22 TRUE 0
## 23 TRUE 0
## 24 TRUE 0
## 25 TRUE 0
## 26 TRUE 0
## 27 TRUE 0
## 28 TRUE 0
## 29 TRUE 0
## 30 TRUE 0
## 31 TRUE 0
## 32 TRUE 0
## 33 TRUE 0
## 34 TRUE 0
## 35 TRUE 0
## 36 TRUE 0
## 37 TRUE 0
## 38 TRUE 0
## 39 TRUE 0
## 40 TRUE 0
## 41 TRUE 0
## 42 TRUE 0
## 43 TRUE 0
## 44 TRUE 0
## 45 TRUE 0
## 46 TRUE 0
## 47 TRUE 0
## 48 TRUE 0
## 49 TRUE 0
## 50 TRUE 0
## 51 TRUE 0
## 52 TRUE 0
mlr <- hf_lags_table(yy ~ mls(yy, 1, 1) + mls(zz, 0:11, 66, nbetaMT),
start = list(zz = c(2, 1, 5, 0)),
from=c(zz=0),to=list(zz=c(5,56)))
##
## Model selection progress:
##
|
| | 0%
|
|= | 2%
|
|== | 4%
|
|==== | 6%
|
|===== | 8%
|
|====== | 10%
|
|======== | 12%
|
|========= | 13%
|
|========== | 15%
|
|=========== | 17%
|
|============ | 19%
|
|============== | 21%
|
|=============== | 23%
|
|================ | 25%
|
|================== | 27%
|
|=================== | 29%
|
|==================== | 31%
|
|===================== | 33%
|
|====================== | 35%
|
|======================== | 37%
|
|========================= | 38%
|
|========================== | 40%
|
|============================ | 42%
|
|============================= | 44%
|
|============================== | 46%
|
|=============================== | 48%
|
|================================ | 50%
|
|================================== | 52%
|
|=================================== | 54%
|
|==================================== | 56%
|
|====================================== | 58%
|
|======================================= | 60%
|
|======================================== | 62%
|
|========================================= | 63%
|
|========================================== | 65%
|
|============================================ | 67%
|
|============================================= | 69%
|
|============================================== | 71%
|
|================================================ | 73%
|
|================================================= | 75%
|
|================================================== | 77%
|
|=================================================== | 79%
|
|==================================================== | 81%
|
|====================================================== | 83%
|
|======================================================= | 85%
|
|======================================================== | 87%
|
|========================================================== | 88%
|
|=========================================================== | 90%
|
|============================================================ | 92%
|
|============================================================= | 94%
|
|============================================================== | 96%
|
|================================================================ | 98%
|
|=================================================================| 100%
# 34 lag is the best
mlr
## model AIC.restricted
## 1 yy ~ mls(yy, 1, 1) + mls(zz, 0:5, 66, nbetaMT) 247.1044
## 2 yy ~ mls(yy, 1, 1) + mls(zz, 0:6, 66, nbetaMT) 245.2506
## 3 yy ~ mls(yy, 1, 1) + mls(zz, 0:7, 66, nbetaMT) 243.3692
## 4 yy ~ mls(yy, 1, 1) + mls(zz, 0:8, 66, nbetaMT) 247.4860
## 5 yy ~ mls(yy, 1, 1) + mls(zz, 0:9, 66, nbetaMT) 248.3736
## 6 yy ~ mls(yy, 1, 1) + mls(zz, 0:10, 66, nbetaMT) 248.7546
## 7 yy ~ mls(yy, 1, 1) + mls(zz, 0:11, 66, nbetaMT) 249.9835
## 8 yy ~ mls(yy, 1, 1) + mls(zz, 0:12, 66, nbetaMT) 250.1077
## 9 yy ~ mls(yy, 1, 1) + mls(zz, 0:13, 66, nbetaMT) 250.1615
## 10 yy ~ mls(yy, 1, 1) + mls(zz, 0:14, 66, nbetaMT) 249.8992
## 11 yy ~ mls(yy, 1, 1) + mls(zz, 0:15, 66, nbetaMT) 250.0948
## 12 yy ~ mls(yy, 1, 1) + mls(zz, 0:16, 66, nbetaMT) 250.1063
## 13 yy ~ mls(yy, 1, 1) + mls(zz, 0:17, 66, nbetaMT) 250.1198
## 14 yy ~ mls(yy, 1, 1) + mls(zz, 0:18, 66, nbetaMT) 249.7835
## 15 yy ~ mls(yy, 1, 1) + mls(zz, 0:19, 66, nbetaMT) 249.9055
## 16 yy ~ mls(yy, 1, 1) + mls(zz, 0:20, 66, nbetaMT) 249.7262
## 17 yy ~ mls(yy, 1, 1) + mls(zz, 0:21, 66, nbetaMT) 246.5597
## 18 yy ~ mls(yy, 1, 1) + mls(zz, 0:22, 66, nbetaMT) 247.7794
## 19 yy ~ mls(yy, 1, 1) + mls(zz, 0:23, 66, nbetaMT) 248.7517
## 20 yy ~ mls(yy, 1, 1) + mls(zz, 0:24, 66, nbetaMT) 249.4000
## 21 yy ~ mls(yy, 1, 1) + mls(zz, 0:25, 66, nbetaMT) 238.4988
## 22 yy ~ mls(yy, 1, 1) + mls(zz, 0:26, 66, nbetaMT) 240.8199
## 23 yy ~ mls(yy, 1, 1) + mls(zz, 0:27, 66, nbetaMT) 239.2180
## 24 yy ~ mls(yy, 1, 1) + mls(zz, 0:28, 66, nbetaMT) 237.7143
## 25 yy ~ mls(yy, 1, 1) + mls(zz, 0:29, 66, nbetaMT) 236.5302
## 26 yy ~ mls(yy, 1, 1) + mls(zz, 0:30, 66, nbetaMT) 236.3929
## 27 yy ~ mls(yy, 1, 1) + mls(zz, 0:31, 66, nbetaMT) 235.3514
## 28 yy ~ mls(yy, 1, 1) + mls(zz, 0:32, 66, nbetaMT) 230.0699
## 29 yy ~ mls(yy, 1, 1) + mls(zz, 0:33, 66, nbetaMT) 233.7186
## 30 yy ~ mls(yy, 1, 1) + mls(zz, 0:34, 66, nbetaMT) 228.9115
## 31 yy ~ mls(yy, 1, 1) + mls(zz, 0:35, 66, nbetaMT) 231.2000
## 32 yy ~ mls(yy, 1, 1) + mls(zz, 0:36, 66, nbetaMT) 231.8221
## 33 yy ~ mls(yy, 1, 1) + mls(zz, 0:37, 66, nbetaMT) 230.1454
## 34 yy ~ mls(yy, 1, 1) + mls(zz, 0:38, 66, nbetaMT) 231.4586
## 35 yy ~ mls(yy, 1, 1) + mls(zz, 0:39, 66, nbetaMT) 236.6759
## 36 yy ~ mls(yy, 1, 1) + mls(zz, 0:40, 66, nbetaMT) 233.5961
## 37 yy ~ mls(yy, 1, 1) + mls(zz, 0:41, 66, nbetaMT) 233.4720
## 38 yy ~ mls(yy, 1, 1) + mls(zz, 0:42, 66, nbetaMT) 232.7385
## 39 yy ~ mls(yy, 1, 1) + mls(zz, 0:43, 66, nbetaMT) 233.7704
## 40 yy ~ mls(yy, 1, 1) + mls(zz, 0:44, 66, nbetaMT) 234.2126
## 41 yy ~ mls(yy, 1, 1) + mls(zz, 0:45, 66, nbetaMT) 234.4101
## 42 yy ~ mls(yy, 1, 1) + mls(zz, 0:46, 66, nbetaMT) 233.1887
## 43 yy ~ mls(yy, 1, 1) + mls(zz, 0:47, 66, nbetaMT) 232.8973
## 44 yy ~ mls(yy, 1, 1) + mls(zz, 0:48, 66, nbetaMT) 234.3591
## 45 yy ~ mls(yy, 1, 1) + mls(zz, 0:49, 66, nbetaMT) 232.8798
## 46 yy ~ mls(yy, 1, 1) + mls(zz, 0:50, 66, nbetaMT) 232.6653
## 47 yy ~ mls(yy, 1, 1) + mls(zz, 0:51, 66, nbetaMT) 237.9361
## 48 yy ~ mls(yy, 1, 1) + mls(zz, 0:52, 66, nbetaMT) 237.7985
## 49 yy ~ mls(yy, 1, 1) + mls(zz, 0:53, 66, nbetaMT) 233.3146
## 50 yy ~ mls(yy, 1, 1) + mls(zz, 0:54, 66, nbetaMT) 232.4974
## 51 yy ~ mls(yy, 1, 1) + mls(zz, 0:55, 66, nbetaMT) 233.0124
## 52 yy ~ mls(yy, 1, 1) + mls(zz, 0:56, 66, nbetaMT) 232.6743
## BIC.restricted AIC.unrestricted BIC.unrestricted hAh.test.p.value First
## 1 266.7326 245.0055 270.2417 0.098249561 FALSE
## 2 264.8787 245.4411 273.4813 0.155533012 FALSE
## 3 262.9973 247.2446 278.0889 0.426520632 FALSE
## 4 267.1141 248.0420 281.6903 0.211649099 FALSE
## 5 268.0017 249.0910 285.5433 0.186005443 FALSE
## 6 268.3828 251.0383 290.2946 0.259689879 FALSE
## 7 269.6117 250.7904 292.8507 0.073724359 FALSE
## 8 269.7358 252.6972 297.5616 0.108582920 FALSE
## 9 269.7897 254.4659 302.1342 0.149217007 FALSE
## 10 269.5273 255.6085 306.0809 0.182503128 FALSE
## 11 269.7229 256.4815 309.7579 0.184348747 FALSE
## 12 269.7345 258.4765 314.5569 0.248409876 FALSE
## 13 269.7479 260.2902 319.1747 0.310270640 FALSE
## 14 269.4116 250.8327 312.5212 0.028330373 FALSE
## 15 269.5337 248.8087 313.3012 0.012818221 FALSE
## 16 269.3543 250.0395 317.3360 0.017434472 FALSE
## 17 266.1878 247.2083 317.3088 0.017088831 FALSE
## 18 267.4075 248.7649 321.6695 0.017871819 FALSE
## 19 268.3799 249.0108 324.7194 0.011812479 FALSE
## 20 269.0282 248.8252 327.3378 0.008279231 FALSE
## 21 258.1270 249.2644 330.5811 0.142145269 FALSE
## 22 260.4481 249.7398 333.8604 0.081260721 FALSE
## 23 258.8462 251.5307 338.4554 0.157411351 FALSE
## 24 257.3425 253.5218 343.2505 0.283779386 FALSE
## 25 256.1584 255.0657 347.5984 0.445372636 FALSE
## 26 256.0210 257.0655 352.4022 0.551970004 FALSE
## 27 254.9795 259.0387 357.1794 0.695326239 FALSE
## 28 249.6981 256.6729 357.6177 0.650314585 FALSE
## 29 253.3468 258.6673 362.4160 0.682300772 FALSE
## 30 248.5397 260.1923 366.7451 0.796535231 FALSE
## 31 250.8282 260.2600 369.6168 0.676506211 FALSE
## 32 251.4502 260.0672 372.2281 0.612086059 FALSE
## 33 249.7736 257.6398 372.6046 0.551634845 FALSE
## 34 251.0868 258.0838 375.8527 0.488160470 FALSE
## 35 256.3040 260.0118 380.5847 0.517555921 FALSE
## 36 253.2242 261.7802 385.1571 0.515519996 FALSE
## 37 253.1002 263.7217 389.9027 0.586165841 FALSE
## 38 252.3666 265.7180 394.7030 0.681074972 FALSE
## 39 253.3986 265.7014 397.4904 0.619612189 FALSE
## 40 253.8407 267.6868 402.2798 0.664939477 FALSE
## 41 254.0383 267.8046 405.2017 0.645404693 FALSE
## 42 252.8168 269.7455 409.9466 0.748399364 FALSE
## 43 252.5255 271.6049 414.6100 0.804465119 FALSE
## 44 253.9872 273.3985 419.2076 0.802491818 FALSE
## 45 252.5079 275.1690 423.7821 0.876250251 FALSE
## 46 252.2935 276.6666 428.0837 0.902831408 FALSE
## 47 257.5643 278.4310 432.6521 0.895795341 FALSE
## 48 257.4267 280.3540 437.3792 0.909498058 FALSE
## 49 252.9427 275.5277 435.3569 0.843590414 FALSE
## 50 252.1256 269.8066 432.4398 0.692619016 FALSE
## 51 252.6406 271.6558 437.0931 0.730920413 FALSE
## 52 252.3025 273.6310 441.8722 0.794763639 FALSE
## Second Convergence
## 1 FALSE 0
## 2 FALSE 0
## 3 FALSE 0
## 4 FALSE 0
## 5 FALSE 0
## 6 FALSE 0
## 7 FALSE 0
## 8 FALSE 0
## 9 FALSE 0
## 10 FALSE 0
## 11 FALSE 0
## 12 FALSE 0
## 13 FALSE 0
## 14 FALSE 0
## 15 FALSE 0
## 16 FALSE 0
## 17 FALSE 0
## 18 FALSE 0
## 19 FALSE 0
## 20 FALSE 0
## 21 TRUE 0
## 22 FALSE 0
## 23 TRUE 0
## 24 FALSE 0
## 25 FALSE 0
## 26 FALSE 0
## 27 FALSE 0
## 28 TRUE 0
## 29 FALSE 0
## 30 TRUE 0
## 31 TRUE 0
## 32 TRUE 0
## 33 TRUE 0
## 34 TRUE 0
## 35 FALSE 0
## 36 TRUE 0
## 37 TRUE 0
## 38 TRUE 0
## 39 TRUE 0
## 40 TRUE 0
## 41 TRUE 0
## 42 TRUE 0
## 43 TRUE 0
## 44 TRUE 0
## 45 TRUE 0
## 46 TRUE 0
## 47 TRUE 0
## 48 FALSE 0
## 49 TRUE 0
## 50 TRUE 0
## 51 TRUE 0
## 52 TRUE 0
mlr <- hf_lags_table(yy~trend+fmls(zz,11,66,nealmon),
start=list(zz=rep(0,3)),
from=c(zz=0),to=list(zz=c(5,56)))
##
## Model selection progress:
##
|
| | 0%
|
|= | 2%
|
|== | 4%
|
|==== | 6%
|
|===== | 8%
|
|====== | 10%
|
|======== | 12%
|
|========= | 13%
|
|========== | 15%
|
|=========== | 17%
|
|============ | 19%
|
|============== | 21%
|
|=============== | 23%
|
|================ | 25%
|
|================== | 27%
|
|=================== | 29%
|
|==================== | 31%
|
|===================== | 33%
|
|====================== | 35%
|
|======================== | 37%
|
|========================= | 38%
|
|========================== | 40%
|
|============================ | 42%
|
|============================= | 44%
|
|============================== | 46%
|
|=============================== | 48%
|
|================================ | 50%
|
|================================== | 52%
|
|=================================== | 54%
|
|==================================== | 56%
|
|====================================== | 58%
|
|======================================= | 60%
|
|======================================== | 62%
|
|========================================= | 63%
|
|========================================== | 65%
|
|============================================ | 67%
|
|============================================= | 69%
|
|============================================== | 71%
|
|================================================ | 73%
|
|================================================= | 75%
|
|================================================== | 77%
|
|=================================================== | 79%
|
|==================================================== | 81%
|
|====================================================== | 83%
|
|======================================================= | 85%
|
|======================================================== | 87%
|
|========================================================== | 88%
|
|=========================================================== | 90%
|
|============================================================ | 92%
|
|============================================================= | 94%
|
|============================================================== | 96%
|
|================================================================ | 98%
|
|=================================================================| 100%
# 21 las is the best
mlr
## model AIC.restricted BIC.restricted
## 1 yy ~ trend + mls(zz, 0:5, 66, nealmon) 290.7060 307.5791
## 2 yy ~ trend + mls(zz, 0:6, 66, nealmon) 289.9306 306.8037
## 3 yy ~ trend + mls(zz, 0:7, 66, nealmon) 290.2153 307.0884
## 4 yy ~ trend + mls(zz, 0:8, 66, nealmon) 290.6376 307.5107
## 5 yy ~ trend + mls(zz, 0:9, 66, nealmon) 290.6973 307.5704
## 6 yy ~ trend + mls(zz, 0:10, 66, nealmon) 290.6999 307.5730
## 7 yy ~ trend + mls(zz, 0:11, 66, nealmon) 289.7547 306.6278
## 8 yy ~ trend + mls(zz, 0:12, 66, nealmon) 289.1709 306.0440
## 9 yy ~ trend + mls(zz, 0:13, 66, nealmon) 289.4602 306.3333
## 10 yy ~ trend + mls(zz, 0:14, 66, nealmon) 289.5182 306.3913
## 11 yy ~ trend + mls(zz, 0:15, 66, nealmon) 289.5178 306.3909
## 12 yy ~ trend + mls(zz, 0:16, 66, nealmon) 289.5182 306.3913
## 13 yy ~ trend + mls(zz, 0:17, 66, nealmon) 289.5178 306.3909
## 14 yy ~ trend + mls(zz, 0:18, 66, nealmon) 288.7308 305.6039
## 15 yy ~ trend + mls(zz, 0:19, 66, nealmon) 289.5177 306.3909
## 16 yy ~ trend + mls(zz, 0:20, 66, nealmon) 288.8028 305.6759
## 17 yy ~ trend + mls(zz, 0:21, 66, nealmon) 276.2289 293.1020
## 18 yy ~ trend + mls(zz, 0:22, 66, nealmon) 282.9951 299.8682
## 19 yy ~ trend + mls(zz, 0:23, 66, nealmon) 278.6477 295.5208
## 20 yy ~ trend + mls(zz, 0:24, 66, nealmon) 278.6408 295.5139
## 21 yy ~ trend + mls(zz, 0:25, 66, nealmon) 278.6571 295.5302
## 22 yy ~ trend + mls(zz, 0:26, 66, nealmon) 278.6448 295.5179
## 23 yy ~ trend + mls(zz, 0:27, 66, nealmon) 278.6519 295.5250
## 24 yy ~ trend + mls(zz, 0:28, 66, nealmon) 288.2292 305.1024
## 25 yy ~ trend + mls(zz, 0:29, 66, nealmon) 287.2526 304.1257
## 26 yy ~ trend + mls(zz, 0:30, 66, nealmon) 287.8846 304.7577
## 27 yy ~ trend + mls(zz, 0:31, 66, nealmon) 285.8887 302.7618
## 28 yy ~ trend + mls(zz, 0:32, 66, nealmon) 282.3834 299.2565
## 29 yy ~ trend + mls(zz, 0:33, 66, nealmon) 284.1042 300.9773
## 30 yy ~ trend + mls(zz, 0:34, 66, nealmon) 277.5367 294.4098
## 31 yy ~ trend + mls(zz, 0:35, 66, nealmon) 281.6342 298.5074
## 32 yy ~ trend + mls(zz, 0:36, 66, nealmon) 282.8750 299.7481
## 33 yy ~ trend + mls(zz, 0:37, 66, nealmon) 282.7329 299.6060
## 34 yy ~ trend + mls(zz, 0:38, 66, nealmon) 282.9941 299.8672
## 35 yy ~ trend + mls(zz, 0:39, 66, nealmon) 287.8524 304.7255
## 36 yy ~ trend + mls(zz, 0:40, 66, nealmon) 284.2225 301.0956
## 37 yy ~ trend + mls(zz, 0:41, 66, nealmon) 284.3645 301.2376
## 38 yy ~ trend + mls(zz, 0:42, 66, nealmon) 284.3166 301.1897
## 39 yy ~ trend + mls(zz, 0:43, 66, nealmon) 284.4349 301.3080
## 40 yy ~ trend + mls(zz, 0:44, 66, nealmon) 284.4421 301.3152
## 41 yy ~ trend + mls(zz, 0:45, 66, nealmon) 284.4246 301.2977
## 42 yy ~ trend + mls(zz, 0:46, 66, nealmon) 284.4273 301.3004
## 43 yy ~ trend + mls(zz, 0:47, 66, nealmon) 284.3988 301.2719
## 44 yy ~ trend + mls(zz, 0:48, 66, nealmon) 284.4420 301.3151
## 45 yy ~ trend + mls(zz, 0:49, 66, nealmon) 284.4889 301.3620
## 46 yy ~ trend + mls(zz, 0:50, 66, nealmon) 284.4873 301.3604
## 47 yy ~ trend + mls(zz, 0:51, 66, nealmon) 284.3665 301.2396
## 48 yy ~ trend + mls(zz, 0:52, 66, nealmon) 284.3699 301.2430
## 49 yy ~ trend + mls(zz, 0:53, 66, nealmon) 284.3744 301.2475
## 50 yy ~ trend + mls(zz, 0:54, 66, nealmon) 284.4813 301.3544
## 51 yy ~ trend + mls(zz, 0:55, 66, nealmon) 284.4093 301.2824
## 52 yy ~ trend + mls(zz, 0:56, 66, nealmon) 284.3613 301.2345
## AIC.unrestricted BIC.unrestricted hAh.test.p.value First Second
## 1 293.4148 318.7244 0.3736772 FALSE TRUE
## 2 294.7929 322.9148 0.5669841 FALSE TRUE
## 3 296.7440 327.6780 0.6638896 FALSE TRUE
## 4 297.2791 331.0253 0.5454724 FALSE TRUE
## 5 296.2509 332.8093 0.3423799 FALSE TRUE
## 6 298.0063 337.3769 0.4280106 FALSE TRUE
## 7 298.1888 340.3716 0.4544803 FALSE TRUE
## 8 299.6978 344.6927 0.5659402 FALSE TRUE
## 9 301.6096 349.4167 0.6291162 FALSE TRUE
## 10 303.6014 354.2208 0.7114354 FALSE TRUE
## 11 305.2633 358.6948 0.7631240 FALSE TRUE
## 12 307.2105 363.4542 0.8252434 FALSE TRUE
## 13 309.1492 368.2050 0.8744073 FALSE TRUE
## 14 305.7268 367.5949 0.6556558 FALSE FALSE
## 15 306.1741 370.8543 0.5742752 FALSE TRUE
## 16 304.0485 371.5409 0.4308590 FALSE TRUE
## 17 301.4316 371.7362 0.9289421 FALSE TRUE
## 18 303.4311 376.5479 0.6633882 FALSE TRUE
## 19 301.8138 377.7428 0.7668234 NA NA
## 20 302.7150 381.4562 0.7698524 NA NA
## 21 304.1010 385.6543 0.7952664 NA NA
## 22 303.6805 388.0460 0.7360691 NA NA
## 23 305.6009 392.7786 0.7889615 NA NA
## 24 307.1192 397.1091 0.3382150 FALSE TRUE
## 25 307.9372 400.7393 0.3951971 FALSE TRUE
## 26 307.0213 402.6355 0.2975731 FALSE TRUE
## 27 306.9483 405.3748 0.3580282 FALSE TRUE
## 28 300.6723 401.9110 0.2242712 FALSE TRUE
## 29 302.4867 406.5375 0.2119644 FALSE TRUE
## 30 301.0772 407.9402 0.4380681 FALSE TRUE
## 31 302.3023 411.9775 0.2564660 FALSE TRUE
## 32 304.2723 416.7596 0.2706747 FALSE TRUE
## 33 303.6634 418.9629 0.2392707 FALSE TRUE
## 34 303.6357 421.7474 0.2164905 FALSE TRUE
## 35 305.5967 426.5207 0.3870314 FALSE TRUE
## 36 307.5800 431.3161 0.2777065 FALSE TRUE
## 37 309.2379 435.7862 0.3178850 FALSE TRUE
## 38 311.2018 440.5623 0.3811224 FALSE TRUE
## 39 310.5370 442.7097 0.3403371 FALSE TRUE
## 40 312.4771 447.4619 0.4004307 FALSE TRUE
## 41 312.6798 450.4768 0.3909288 FALSE TRUE
## 42 313.7726 454.3818 0.4263544 FALSE TRUE
## 43 315.1721 458.5935 0.4711770 FALSE TRUE
## 44 314.9725 461.2061 0.4487464 FALSE TRUE
## 45 316.9097 465.9554 0.5210889 FALSE TRUE
## 46 318.8954 470.7534 0.5886732 FALSE TRUE
## 47 320.8913 475.5615 0.6467043 FALSE TRUE
## 48 322.0095 479.4918 0.6778325 FALSE TRUE
## 49 322.4886 482.7831 0.6873166 FALSE TRUE
## 50 319.7858 482.8925 0.5892686 FALSE TRUE
## 51 321.6549 487.5738 0.6475066 FALSE TRUE
## 52 323.6420 492.3730 0.7075794 FALSE TRUE
## Convergence
## 1 0
## 2 0
## 3 0
## 4 0
## 5 0
## 6 0
## 7 0
## 8 0
## 9 0
## 10 0
## 11 0
## 12 0
## 13 0
## 14 0
## 15 0
## 16 0
## 17 0
## 18 0
## 19 0
## 20 0
## 21 1
## 22 0
## 23 1
## 24 0
## 25 0
## 26 0
## 27 0
## 28 0
## 29 0
## 30 0
## 31 0
## 32 0
## 33 0
## 34 0
## 35 0
## 36 0
## 37 0
## 38 0
## 39 0
## 40 0
## 41 0
## 42 0
## 43 0
## 44 0
## 45 0
## 46 0
## 47 0
## 48 0
## 49 0
## 50 0
## 51 0
## 52 0
Use lapply function to test how many lags are the best choice.
#
#
# trend <- 1:length(yy)
# #allk <- lapply(c(12,15,18,24)-1, function(k) {
# #midas_r(midas_r(yy~trend+fmls(xx,k,3,nealmon),start=list(xx=rep(0,3))), Ofunction="nls")
# #})
#
# # test K for xx
#
# allk <- lapply(c(3,6,9,12)-1, function(k) {
# midas_r(midas_r(yy~trend+fmls(xx,k,3,nealmon),
# start=list(xx=rep(0,3))), Ofunction="nls")
# })
#
#
# ####Compute the derivative test
# dtest <- lapply(allk,deriv_tests)
#
#
# ###The first derivative tests, gradient is zero
# sapply(dtest,with,first)
#
# ###The second derivative tests, hessian is positive definite
# sapply(dtest,with,second)
#
#
# ###The minimal eigenvalue of hessian is borderline zero, yet positive.
# sapply(dtest,with,min(eigenval))
#
#
# ###Apply hAh test
# lapply(allk,hAh.test)
#
# ###Apply robust hAh test
# lapply(allk,hAhr.test)
# ###View summaries
# lapply(allk,summary)
# lapply(allk,coef)
#
#
# ##Plot the coefficients
# #dev.new()
# par(mfrow=c(2,2))
#
# lapply(allk,function(x){
# cfur <- coef(x$unrestricted)
# cfur <- cfur[grep("fmls",names(cfur))]
# cfre <- weight_coef(x)
# k <- length(cfur)
# sdval <- sqrt(diag(vcovHAC(x$unrestricted)))
# sdval <- sdval[grep("fmls",names(sdval))]
# plot(0:(k-1),cfur,col="black",ylab="Beta coefficients",xlab="h")
# title(main=sprintf("d = %.0f: p-val.(hAh_HAC) < %.2f", k,
# max(hAhr.test(x)$p.value, 0.01)),
# cex.main = 1, font.main = 4, col.main = "black")
# points(c(0:(k - 1)), cfre, type = "l", col = "blue")
# points(c(0:(k - 1)), cfur + 2 * sdval[1:k], type = "l", col = "red", lty = 2)
# points(c(0:(k - 1)), cfur - 2 * sdval[1:k], type = "l", col = "red", lty = 2)
# })
#
#https://github.com/mpiktas/midasr/blob/d8bae1f1f962f71da9cee517415623f721a7a8cf/R/lagspec.R
##' normalized exponential Almon lag restricts the coefficients \eqn{theta_h} in the following way:
##'
##' \deqn{\theta_{h}=\delta\frac{\exp(\lambda_1(h+1)+\dots+\lambda_r(h+1)^r)}{\sum_{s=0}^d\exp(\lambda_1(s+1)+\dots+\lambda_r(h+1)^r)}}
##'
##' The parameter \eqn{\delta} should be the first element in vector \code{p}. The degree of the polynomial is then decided by the number of the remaining parameters.
##' @export
nealmon <- function(p,d,m) {
i <- 1:d
plc <- poly(i,degree=length(p)-1,raw=TRUE) %*% p[-1]
as.vector(p[1] * exp(plc)/sum(exp(plc)))
}
##' Normalized beta probability density function MIDAS weights specification (MATLAB toolbox compatible)
##' Calculate MIDAS weights according to normalized beta probability density function specification. Compatible with the specification in MATLAB toolbox.
##' @param p parameters for normalized beta probability density function
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @export
nbetaMT <- function(p,d,m) {
eps <- .Machine$double.eps
xi <- (1:d - 1)/(d - 1)
xi[1] <- xi[1]+eps
xi[d] <- xi[d]-eps
nb <- xi^(p[2] - 1) * (1 - xi)^(p[3] - 1)
if(sum(nb)<eps) {
if(abs(p[4])<eps) rep(0,d)
else p[1]*rep(1/d,length(nb))
} else {
w <- (nb/sum(nb) + p[4])
p[1]*w/sum(w)
}
}
##' Normalized beta probability density function MIDAS weights specification
##' Calculate MIDAS weights according to normalized beta probability density function specification
##' @param p parameters for normalized beta probability density function
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @export
nbeta <- function(p,d,m) {
nbetaMT(c(p,0),d,m)
}