# R/lagspec.R In midasr: Mixed Data Sampling Regression

##' Normalized Exponential Almon lag MIDAS coefficients
##'
##' Calculate normalized exponential Almon lag coefficients given the parameters and required number of coefficients.
##' @param p parameters for Almon lag
##' @param d number of the coefficients
##' @param m the frequency, currently ignored.
##' @return vector of coefficients
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @examples
##'
##' data("USunempr")
##' data("USrealgdp")
##'
##' y <- diff(log(USrealgdp))
##' x <- window(diff(USunempr),start=1949)
##' t <- 1:length(y)
##'
##' midas_r(y~t+fmls(x,11,12,nealmon),start=list(x=c(0,0,0)))
##'
##' @details Given unrestricted MIDAS regression
##'
##' \deqn{y_t=\sum_{h=0}^d\theta_{h}x_{tm-h}+\mathbf{z_t}\beta+u_t}
##'
##' 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.
##' @importFrom stats poly
##' @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)))
}

##' Produces weights for aggregates based MIDAS regression
##'
##' Suppose a weight function \eqn{w(\beta,\theta)} satisfies the following equation:
##' \deqn{w(\beta,\theta)=\beta g(\theta)}
##'
##' The following combinations are defined, corresponding to structure types \code{A}, \code{B} and \code{C} respectively:
##' \deqn{(w(\beta_1,\theta_1),...,w(\beta_k,\theta_k))}
##' \deqn{(w(\beta_1,\theta),...,w(\beta_k,\theta))}
##' \deqn{\beta(w(1,\theta),...,w(1,\theta)),}
##'
##' where \eqn{k} is the number of low frequency lags, i.e. \eqn{d/m}. If the latter
##' value is not whole number, the error is produced.
##'
##'
##' The starting values \eqn{p} should be supplied then as follows:
##' \deqn{(\beta_1,\theta_1,...,\beta_k,\theta_k)}
##' \deqn{(\beta_1,...,\beta_k,\theta)}
##' \deqn{(\beta,\theta)}
##'
##'
##' @title Weights for aggregates based MIDAS regressions
##' @param p parameters for weight functions, see details.
##' @param d number of high frequency lags
##' @param m the frequency
##' @param weight the weight function
##' @param type type of structure, a string, one of A, B or C.
##' @return a vector of weights
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @export
amweights <- function(p,d,m,weight=nealmon,type=c("A","B","C")) {
type <- match.arg(type)
hf <- d%/%m

if(d%%m!=0)stop("Number of high frequency lags should be a multiple of frequency")
if(type=="A") {
return(as.vector(apply(matrix(1:length(p),ncol=hf),2,function(x)weight(p[x],m,m))))
}
if(type=="B") {
theta <- p[-1:-hf]
return(as.vector(sapply(p[1:hf],function(beta)weight(c(beta,theta),m))))
}
if(type=="C") {
return(p[1]*as.vector(apply(matrix(2:length(p),ncol=hf),2,function(x)weight(c(1,p[x]),m,m))))
}
}

##' Gradient function for normalized exponential Almon lag weights
##'
##' Gradient function for normalized exponential Almon lag weights
##' @param p hyperparameters for Almon lag
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @author Vaidotas Zemlys
##' @export
i <- 1:d
pl <- poly(i,degree=length(p)-1,raw=TRUE)
eplc <- exp(pl%*%p[-1])[,,drop=TRUE]
ds <- apply(pl*eplc,2,sum)
s <- sum(eplc)
cbind(eplc/s,p[1]*(pl*eplc/s-eplc%*%t(ds)/s^2))
}

##' 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[1:3],0),d,m)
}

##' Gradient function for normalized beta probability density function MIDAS weights specification
##'
##' Calculate gradient function for normalized beta probability density function specification of MIDAS weights.
##' @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
}

##' 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) } } ##' Gradient function for normalized beta probability density function MIDAS weights specification (MATLAB toolbox compatible) ##' ##' Calculate gradient function for normalized beta probability density function specification of MIDAS weights. ##' @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_gradient <- 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)
snb <- sum(nb)
wnb <- 1+d*p[4]
if(snb>eps) {
nba <- nb*log(xi)
nbb <- nb*log(1-xi)
a <- nba/snb-nb*sum(nba)/snb^2
b <- nbb/snb-nb*sum(nbb)/snb^2
cbind((nb/snb+p[4])/wnb,p[1]*a/wnb,p[1]*b/wnb,p[1]*(1-nb/snb*d)/wnb^2)
}
else {
gres <- matrix(0,nrow=d,ncol=4)
if(abs(p[4])>eps)gres[,1] <- 1/d
gres
}
}

##' Almon polynomial MIDAS weights specification
##'
##' Calculate Almon polynomial MIDAS weights
##' @param p parameters for Almon polynomial weights
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @export
almonp <- function(p,d,m) {
i <- 1:d
plc <- poly(i,degree=length(p)-1,raw=TRUE) %*%p[-1]+p[1]
as.vector(plc)
}

##' Gradient function for Almon polynomial MIDAS weights
##'
##' Calculate gradient for Almon polynomial MIDAS weights specification
##' @param p vector of parameters for Almon polynomial specification
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Vaidotas Zemlys
##' @export
i <- 1:d
plc <- poly(i,degree=length(p)-1,raw=TRUE)
cbind(1,plc)
}

##' Step function specification for MIDAS weights
##'
##' Step function specification for MIDAS weights
##' @param p vector of parameters
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @param a vector of increasing positive integers indicating the steps
##' @return vector of coefficients
##' @author Vaidotas Zemlys
##' @export
polystep <- function(p,d,m,a) {
if(length(a)!=length(p)-1)stop("The number of steps should be number of parameters minus one")
if(min(a)<=1 | max(a)>=d)stop("The steps are out of bounds")
a <- c(0,a,d)
rep(p,times=diff(a))
}
##' Gradient of step function specification for MIDAS weights
##'
##' Gradient of step function specification for MIDAS weights
##' @param p vector of parameters
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @param a vector of increasing positive integers indicating the steps
##' @return vector of coefficients
##' @author Vaidotas Zemlys
##' @export
if(length(a)!=length(p)-1)stop("The number of steps should be number of parameters minus one")
if(min(a)<=1 | max(a)>=d)stop("The steps are out of bounds")
a <- c(0,a,d)
da <- diff(a)
r <- cbind(a[-length(a)],da)
apply(r,1,function(x) {
res <- rep(0,d)
res[(x[1]+1):(x[1]+x[2])] <- 1
res
})
}

##' Normalized Gompertz probability density function MIDAS weights specification
##'
##' Calculate MIDAS weights according to normalized Gompertz probability density function specification
##' @param p parameters for normalized Gompertz probability density function
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Julius Vainora
##' @export
gompertzp <- function(p, d, m) {
i <- 1:d / d
gm <- exp(p[3] * i - p[2] * exp(p[3] * i))
p[1] * gm / sum(gm)
}

##' Gradient function for normalized Gompertz probability density function MIDAS weights specification
##'
##' Calculate gradient function for normalized Gompertz probability density function specification of MIDAS weights.
##' @param p parameters for normalized Gompertz probability density function
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Julius Vainora
##' @export
gompertzp_gradient <- function(p, d, m) {
i <- 1:d / d
gm <- exp(p[3] * i - p[2] * exp(p[3] * i))
dp2 <- -gm * exp(i * p[3])
dp3 <- -gm * i * (p[2] * exp(i * p[3]) - 1)
cbind(gm, p[1] * (dp2 - gm * sum(dp2) / sum(gm)),
p[1] * (dp3 - gm * sum(dp3) / sum(gm))) / sum(gm)
}

##' Normalized Nakagami probability density function MIDAS weights specification
##'
##' Calculate MIDAS weights according to normalized Nakagami probability density function specification
##' @param p parameters for normalized Nakagami probability density function
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Julius Vainora
##' @export
nakagamip <- function(p, d, m) {
i <- 1:d / d
ng <- i^(2 * p[2] - 1) * exp(-p[2] / p[3] * i^2)
p[1] * ng / sum(ng)
}

##' Gradient function for normalized Nakagami probability density function MIDAS weights specification
##'
##' Calculate gradient function for normalized Nakagami probability density function specification of MIDAS weights.
##' @param p parameters for normalized Nakagami probability density function
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Julius Vainora
##' @export
nakagamip_gradient <- function(p, d, m) {
i <- 1:d / d
ng <- i^(2 * p[2] - 1) * exp(-p[2] / p[3] * i^2)
dp2 <- ((2 * log(i) * p[3] - i^2) / p[3]) * ng
dp3 <- (p[2] * i^2 / p[3]^2) * ng
cbind(ng, (dp2 - ng * sum(dp2) / sum(ng)) * p[1],
(dp3 - ng * sum(dp3) / sum(ng)) * p[1]) / sum(ng)
}

##' Normalized log-Cauchy probability density function MIDAS weights specification
##'
##' Calculate MIDAS weights according to normalized log-Cauchy probability density function specification
##' @param p parameters for normalized log-Cauchy probability density function
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Julius Vainora
##' @export
lcauchyp <- function(p, d, m) {
i <- 1:d / d
lc <- 1 / (i * ((log(i) - p[2])^2 + p[3]^2))
p[1] * lc / sum(lc)
}

##' Gradient function for normalized log-Cauchy probability density function MIDAS weights specification
##'
##' Calculate gradient function for normalized log-Cauchy probability density function specification of MIDAS weights.
##' @param p parameters for normalized log-Cauchy probability density function
##' @param d number of coefficients
##' @param m the frequency ratio, currently ignored
##' @return vector of coefficients
##' @author Julius Vainora
##' @export
lcauchyp_gradient <- function(p, d, m) {
i <- 1:d / d
lc <- 1 / (i * ((log(i) - p[2])^2 + p[3]^2))
dp2 <- 2 * lc^2 * i * (log(i) - p[2])
dp3 <- -2 * lc^2 * p[3] * i
cbind(lc, p[1] * (dp2 - lc * sum(dp2) / sum(lc)),
p[1] * (dp3 - lc * sum(dp3) / sum(lc))) / sum(lc)
}

##' HAR(3)-RV model MIDAS weights specification
##'
##' MIDAS weights for Heterogeneous Autoregressive model of Realized Volatilty (HAR-RV). It is assumed that month has 20 days.
##' @title HAR(3)-RV model MIDAS weights specification
##' @param p parameters for Almon lag
##' @param d number of the coefficients
##' @param m the frequency, currently ignored.
##' @return vector of coefficients
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @references Corsi, F., \emph{A Simple Approximate Long-Memory Model of Realized Volatility}, Journal of Financial Econometrics Vol. 7 No. 2 (2009) 174-196
##' @export
harstep <- function(p,d,m) {
if(d!=20) stop("HAR(3)-RV process requires 20 lags")
out <- rep(0,20)
out[1] <- p[1]+p[2]/5+p[3]/20
out[2:5] <- p[2]/5+p[3]/20
out[6:20] <- p[3]/20
out
}
##' Gradient function for HAR(3)-RV model MIDAS weights specification
##'
##' MIDAS weights for Heterogeneous Autoregressive model of Realized Volatilty (HAR-RV). It is assumed that month has 20 days.
##' @title Gradient function for HAR(3)-RV model MIDAS weights specification
##' @param p parameters for Almon lag
##' @param d number of the coefficients
##' @param m the frequency, currently ignored.
##' @return vector of coefficients
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @references Corsi, F., \emph{A Simple Approximate Long-Memory Model of Realized Volatility}, Journal of Financial Econometrics Vol. 7 No. 2 (2009) 174-196
##' @export
if(d!=20) stop("HAR(3)-RV process requires 20 lags")
out <- matrix(0,ncol=3,nrow=d)
out[1,1] <- 1
out[1:5,2] <- 1/5
out[1:20,3] <- 1/20
out
}

#' Calculates the MIDAS coefficients for generalized exponential MIDAS lag specification
#'
#' Generalized exponential MIDAS lag specification is a generalization of exponential
#' Almon lag. It is defined as a product of first order polynomial with exponent
#' of the second order polynomial. This spefication was used by V. Kvedaras and
#' V. Zemlys (2012).
#'
#' @title Generalized exponential MIDAS coefficients
#' @param p a vector of parameters
#' @param d number of coefficients
#' @param m the frequency, currently ignored
#'
#' @return a vector of coefficients
#' @export
#' @author Virmantas Kvedaras, Vaidotas Zemlys
#' @references Kvedaras V., Zemlys, V. \emph{Testing the functional constraints on parameters in regressions with variables of different frequency} Economics Letters 116 (2012) 250-254
genexp <- function(p, d, m) {
i <- (1:d-1)/100
pol <- p[3]*i + p[4]*i^2
(p[1] + p[2]*i)*exp(pol)
}

#' Calculates the gradient of generalized exponential MIDAS lag specification
#'
#' Generalized exponential MIDAS lag specification is a generalization of exponential
#' Almon lag. It is defined as a product of first order polynomial with exponent
#' of the second order polynomial. This spefication was used by V. Kvedaras and
#' V. Zemlys (2012).
#'
#' @title Gradient of generalized exponential MIDAS coefficient generating function
#' @param p a vector of parameters
#' @param d number of coefficients
#' @param m the frequency, currently ignored
#'
#' @return a vector of coefficients
#' @export
#' @author Virmantas Kvedaras, Vaidotas Zemlys
#' @references Kvedaras V., Zemlys, V. \emph{Testing the functional constraints on parameters in regressions with variables of different frequency} Economics Letters 116 (2012) 250-254
genexp_gradient <- function(p, d, m) {
i <- (1:d-1)/100
pol <- p[3]*i + p[4]*i^2
pol0 <- p[1] + p[2]*i
epl <- exp(pol)
cbind(epl, epl*i, pol0*epl*i, pol0*epl*i^2)
}

is_weight_normalized <- function(wf, init, chk = rnorm(5), tol = 1e-8, nc = 1) {

tst <- lapply(c(0,rnorm(chk)), function(x) {
ss <- init
ss[nc] <- ss[nc] + x
try(abs(sum(wf(ss)) - ss[nc]))
})

fail <- sapply(tst, inherits, "try-error")
got_na <- fail
got_na[!fail] <- sapply(tst[!fail], is.na)

if (all(fail) | all(got_na)) stop("Error in checking normalization, all candidate starting values produced errors")

eq <- unlist(tst[!got_na])
if (all(eq < tol)) return(TRUE)
else return(FALSE)
}

find_normalisation_coefficient <- function(wf, init) {
tst <- lapply(1:length(init), function(i) try(is_weight_normalized(wf, init, nc = i)))
tst <- sapply(tst, function(l) if (inherits(l, "try-error")) return(FALSE) else return(l))
which(tst)
}


## Try the midasr package in your browser

Any scripts or data that you put into this service are public.

midasr documentation built on Feb. 23, 2021, 5:11 p.m.