Nothing
#' Mixed causal-noncausal autoregressions with exogenous regressors.
#'
#' @title The MARX function
#' @description This interface-based function allows you to perform model selection for MARX models based on information criteria.
#' @param y Data vector of time series observations.
#' @param x Matrix of data (every column represents one time series). Specify NULL or "not" if not wanted.
#' @param p_max Maximum number of autoregressive parameters (leads + lags) to be included.
#' @param sig_level Significance level for the construction of inference.
#' @param p_C Number of lags (if not specified by the user a model selection procedure is used to determine the number of lags).
#' @param p_NC Number of leads (if not specified by the user a model selection procedure is used to determine the number of leads).
#' @keywords estimation
#' @keywords selection
#' @return The function returns the values of the information criteria for the pseudo-causal models. The user is asked to choose a value for "p". Extensive output for the MARX(r,s,q) model (with p = r + s) which maximizes the log-likelihood is reported.
#' @author Sean Telg
#' @export
#' @examples
#' data <- sim.marx(c('t',1,1), c('t',1,1),100,0.5,0.4,0.3)
#' p_max <- 8
#' sig_level <- 0.05
#' marx(data$y, data$x, p_max, sig_level,1,1) ## p_C and p_NC chosen to be 1: MARX(1,1,1) output.
#' marx(data$y, NULL, p_max,sig_level,1,1) ## MAR(1,1), no exogenous variable specified.
marx <- function(y,x,p_max,sig_level,p_C,p_NC){
if (is.null(x)){
x <- "not"
}
print(match.call())
nargin <- length(as.list(match.call())) -1
if (nargin == 4){
p = 1
}
else{
p = 0
}
if (p == 1){
selection.lag(y,x,p_max)
p_pseudo <- readline(prompt = "Choose lag order for pseudo causal model: ")
p_pseudo <- as.numeric(p_pseudo)
pseudo <- arx.ls(y,x,p_pseudo)
Cov_pseudo <- pseudo[[4]]
U_pseudo <- pseudo[[5]]
test_cdf_pseudo <- cbind(U_pseudo, stats::pnorm(U_pseudo,0,Cov_pseudo))
kstest_results <- stats::ks.test(test_cdf_pseudo[,1],"pnorm",0,Cov_pseudo)
jarquebera <- tseries::jarque.bera.test(U_pseudo)
if (kstest_results$p.value < 0.05){
hh_pseudo = 1 ## reject
}
else{
hh_pseudo = 0 ## not reject
}
if (jarquebera$p.value < 0.05){
jarque_check = 1
}
else{
jarque_check = 0
}
if (hh_pseudo == 0){
cat(' ', "\n")
cat(' ', "\n")
cat('THE KS-TEST FAILS TO REJECT THE NULL OF NORMALLY DISTRIBUTED RESIDUALS OF THE PURELY CAUSAL ARX MODEL', "\n")
cat('p-value:')
cat(kstest_results$p.value, "\n")
cat('WARNING: MIxED ARX MODEL MIGHT NOT BE IDENTIFIABLE!', "\n")
cat(' ', "\n")
cat(' ', "\n")
}
else{
cat(' ', "\n")
cat(' ', "\n")
cat('THE KS-TEST REJECTS THE NULL OF NORMALLY DISTRIBUTED RESIDUALS OF THE PURELY CAUSAL ARX MODEL', "\n")
cat('p-value:')
cat(kstest_results$p.value, "\n")
cat(' ', "\n")
cat(' ', "\n")
}
if (jarque_check == 0){
cat(' ', "\n")
cat(' ', "\n")
cat('THE JB-TEST FAILS TO REJECT THE NULL OF NORMALLY DISTRIBUTED RESIDUALS OF THE PURELY CAUSAL ARX MODEL', "\n")
cat('p-value:')
cat(jarquebera$p.value, "\n")
cat('WARNING: MIxED ARX MODEL MIGHT NOT BE IDENTIFIABLE!', "\n")
cat(' ', "\n")
cat(' ', "\n")
}
else{
cat(' ', "\n")
cat(' ', "\n")
cat('THE JB-TEST REJECTS THE NULL OF NORMALLY DISTRIBUTED RESIDUALS OF THE PURELY CAUSAL ARX MODEL', "\n")
cat('p-value:')
cat(jarquebera$p.value, "\n")
cat(' ', "\n")
cat(' ', "\n")
}
stats::qqnorm(U_pseudo, main="Normal Probability Plot of Residuals")
stats::qqline(U_pseudo)
selection.lag.lead_results <- selection.lag.lead(y,x,p_pseudo)
p_C <- selection.lag.lead_results[[1]]
p_NC <- selection.lag.lead_results[[2]]
}
marx.t_results <- marx.t(y,x,p_C,p_NC);
B_C <- marx.t_results[[1]]
B_NC <- marx.t_results[[2]]
B_x <- marx.t_results[[3]]
IC <- marx.t_results[[4]]
sig <- marx.t_results[[5]]
df <- marx.t_results[[6]]
inference_results <- inference(y,x,B_C,B_NC,B_x,IC,sig,df,sig_level)
BC <- inference_results[[1]]
BNC <- inference_results[[2]]
Bx <- inference_results[[3]]
ICC <- inference_results[[4]]
std_c <- inference_results[[5]]
std_nc <- inference_results[[6]]
std_x <- inference_results[[7]]
std_ic <- inference_results[[8]]
cat(' ', "\n")
cat(' ', "\n")
cat('Estimated Causal Parameters: ', "\n")
print(B_C)
cat(' ', "\n")
cat('Corresponding Standard Errors: ', "\n")
print(std_c)
cat(' ', "\n")
cat('Corresponding Confidence Intervals: ', "\n")
print(BC)
cat(' ', "\n")
cat(' ', "\n")
cat('Estimated Noncausal Parameters: ', "\n")
print(B_NC)
cat(' ', "\n")
cat('Corresponding Standard Errors: ', "\n")
print(std_nc)
cat(' ', "\n")
cat('Corresponding Confidence Intervals: ', "\n")
print(BNC)
cat(' ', "\n")
cat(' ', "\n")
cat('Estimated Parameters Exogenous Variables: ', "\n")
print(B_x)
cat(' ', "\n")
cat('Corresponding Standard Errors: ', "\n")
print(std_x)
cat(' ', "\n")
cat('Corresponding Confidence Intervals: ', "\n")
print(Bx)
cat(' ', "\n")
cat(' ', "\n")
cat('Estimated Intercept: ', "\n")
print(IC)
cat(' ', "\n")
cat('Corresponding Standard Errors: ', "\n")
print(std_ic)
cat(' ', "\n")
cat('Corresponding Confidence Intervals: ', "\n")
print(ICC)
cat(' ', "\n")
cat(' ', "\n")
cat('Estimated Distributional Parameters (df,sig): ', "\n")
print(c(df,sig))
cat(' ', "\n")
cat(' ', "\n")
}
#' @title The regressor matrix function
#' @description This function allows you to create a regressor matrix.
#' @param y Data vector of time series observations.
#' @param x Matrix of data (every column represents one time series). Specify NULL or "not" if not wanted.
#' @param p Number of autoregressive terms to be included.
#' @keywords estimation
#' @return \item{Z}{Regressor matrix}
#' @author Sean Telg
#' @export
#' @examples
#' data <- sim.marx(c('t',3,1),c('t',1,1),100,0.5,0.4,0.3)
#' regressor.matrix(data$y, data$x, 2)
regressor.matrix <- function(y,x,p){
if (is.null(x)){
x <- "not"
}
y <- fBasics::vec(y)
n <- length(y)
if (p==1){
k <-1
}
else{
k <- NCOL(y)
}
if (p > 0){
Z <- matlab::zeros(n,k*p)
for (i in 1:p){
Z[(1+i):n,((i-1)*k+1):(i*k)] <- y[1:(n-i)]
}
Z <- Z[(1+p):n,]
}
else{
Z <- matrix(,nrow=n,ncol=0)
}
if (x == "not" && length(x) == 1){
Z <- Z
}
if (NCOL(x) == 1 && x != "not"){
Z <- cbind(Z,x[(1+p):n])
}
else if (NCOL(x) > 1 && x != "not"){
Z <- cbind(Z,x[(1+p):n,])
}
return(matrix = Z)
}
#' @title The ARX estimation by OLS function
#' @description This function allows you to estimate ARX models by ordinary least squares (OLS).
#' @param y Data vector of time series observations.
#' @param x Matrix of data (every column represents one time series). Specify NULL or "not" if not wanted.
#' @param p Number of autoregressive terms to be included.
#' @keywords estimation
#' @keywords pseudo-causal
#' @return \item{coefficients}{Vector of estimated coefficients.}
#' @return \item{coef.auto}{Vector of estimated autoregressive parameters.}
#' @return \item{coef.exo}{Vector of estimated exogenous parameters.}
#' @return \item{mse}{Mean squared error.}
#' @return \item{residuals}{Residuals.}
#' @return \item{loglikelihood}{Value of the loglikelihood.}
#' @return \item{fitted.values}{Fitted values.}
#' @return \item{df}{Degrees of freedom.}
#' @return \item{vcov}{Variance-covariance matrix of residuals.}
#' @author Sean Telg
#' @export
#' @examples
#' data <- sim.marx(c('t',3,1),c('t',1,1),100,0.5,0.4,0.3)
#' arx.ls(data$y,data$x,2)
arx.ls <- function(y,x,p){
if (is.null(x)){
x <- "not"
}
n <- length(y) - p
Y <- y[(p+1):length(y)]
int <- rep(1,(length(y)-p))
Z <- regressor.matrix(y,x,p)
Z <- cbind(int,Z)
df <- nrow(Z) - NCOL(Z)
B <- solve(t(Z) %*% Z) %*% (t(Z) %*% Y)
if (p > 0){
if (length(x) > 1){
rownames(B) <- c('int', paste('lag', 1:p), paste('exo', 1:NCOL(x)))
}
else{
rownames(B) <- c('int', paste('lag', 1:p))
}
}
else{
if (length(x) > 1){
rownames(B) <- c('int', paste('exo', 1:NCOL(x)))
}
else{
rownames(B) <- 'int'
}
}
FV <- Z %*% B
U <- Y - FV
sig <- (t(U) %*% U)
sig <- as.numeric(sig)
Cov <- (1/n)*sig
Cov <- as.numeric(Cov)
sigma2 <- sum((Y - Z %*% B)^2)/df
qz <- qr(Z)
vcov <- sigma2*chol2inv(qz$qr)
colnames(vcov) <- rownames(vcov) <- colnames(Z)
Loglik <- -(n/2)*(1 + log(2*pi)+log(Cov))
if (p == 0){
B_auto <- 0
}
else{
B_auto <- B[2:(p+1)]
}
if (length(x) > 1){
B_x <- B[(p+2):length(B)]
}
else{
B_x <- 0
}
return(list(coefficients = B, coef.auto = B_auto, coef.exo = B_x, mse = Cov, residuals = U, loglikelihood = Loglik, fitted.values = FV, df = df,vcov=vcov))
}
#' @title The estimation of the MARX model by t-MLE function
#' @description This function allows you to estimate the MARX model by t-MLE.
#' @param y Data vector of time series observations.
#' @param x Matrix of data (every column represents one time series). Specify NULL or "not" if not wanted.
#' @param p_C Number of lags.
#' @param p_NC Number of leads.
#' @param params0 Starting values for the parameters to be estimated (both model and distributional parameters).
#' @keywords estimation
#' @keywords causal-noncausal
#' @return \item{coef.c}{Estimated causal coefficients.}
#' @return \item{coef.nc}{Estimated noncausal coefficients.}
#' @return \item{coef.exo}{Estimated exogenous coefficients.}
#' @return \item{coef.int}{Estimated intercept.}
#' @return \item{scale}{Estimated scale parameter.}
#' @return \item{df}{Estimated degrees of freedom.}
#' @return \item{residuals}{Residuals.}
#' @return \item{se.dist}{Standard errors of the distributional parameters.}
#' @author Sean Telg
#' @export
#' @examples
#' data <- sim.marx(c('t',3,1),c('t',3,1),100,0.5,0.4,0.3)
#' marx.t(data$y,data$x,1,1)
marx.t <- function(y,x,p_C,p_NC,params0){
#print(match.call())
nargin <- length(as.list(match.call())) - 1
if (is.null(x)){
x <- "not"
}
if (length(x) == 1){
numcol <- 0
}
else{
numcol <- NCOL(x)
}
if(numcol > 1){
x.rev <- matrix(data=NA,nrow=length(x[,1]),ncol=numcol)
for (i in 1:numcol){
x.rev[,i] <- rev(x[,i])
}
}
else{
x.rev <- matrix(data=NA,nrow=length(x),ncol=numcol)
x.rev <- rev(x)
}
if (nargin < 5){
y <- fBasics::vec(y)
z <- rev(y)
z <- fBasics::vec(z)
BC0 <- arx.ls(y,x,p_C)[[2]]
Bx0 <- arx.ls(y,x,p_C)[[3]]
BNC0 <- arx.ls(z,x.rev,p_NC)[[2]]
IC0 <- 0
df0 <- 20
sig0 <- 2
BC0 <- fBasics::vec(BC0)
BNC0 <- fBasics::vec(BNC0)
Bx0 <- fBasics::vec(Bx0)
if (length(x) > 1){
if (p_C > 0 && p_NC > 0){
params0 <- rbind(BC0,BNC0,Bx0,IC0,sig0,df0)
}
else if (p_NC > 0 && p_C == 0){
params0 <- rbind(BNC0,Bx0,IC0,sig0,df0)
}
else if (p_C > 0 && p_NC == 0){
params0 <- rbind(BC0,Bx0,IC0,sig0,df0)
}
else if (p_C == 0 && p_NC == 0){
params0 <- rbind(Bx0,IC0,sig0,df0)
}
}
else{
if (p_C > 0 && p_NC > 0){
params0 <- rbind(BC0,BNC0,IC0,sig0,df0)
}
else if (p_NC > 0 && p_C == 0){
params0 <- rbind(BNC0,IC0,sig0,df0)
}
else if (p_C > 0 && p_NC == 0){
params0 <- rbind(BC0,IC0,sig0,df0)
}
else if (p_C == 0 && p_NC == 0){
params0 <- rbind(IC0,sig0,df0)
}
}
}
optimization_results <- stats::optim(params0,ll.max,gr=NULL,y=fBasics::vec(y),p_C=p_C,p_NC=p_NC,x=x,method="BFGS",hessian=TRUE)
PARAMS <- optimization_results$par
if (length(x) > 1){
numcol <- NCOL(x)
if (p_C > 0 && p_NC > 0){
B_C <- PARAMS[1:p_C]
B_NC <- PARAMS[(p_C+1):(p_C + p_NC)]
B_x <- PARAMS[(p_C + p_NC + 1):(p_C + p_NC + numcol)]
IC <- PARAMS[(p_C + p_NC + numcol + 1)]
sig <- PARAMS[(p_C + p_NC + numcol + 2)]
df <- PARAMS[(p_C + p_NC + numcol + 3)]
}
else if (p_NC > 0 && p_C == 0){
B_C <- 0
B_NC <- PARAMS[1:p_NC]
B_x <- PARAMS[(p_NC + 1):(p_NC + numcol)]
IC <- PARAMS[(p_NC + numcol + 1)]
sig <- PARAMS[(p_NC + numcol + 2)]
df <- PARAMS[(p_NC + numcol + 3)]
}
else if (p_C > 0 && p_NC == 0){
B_NC <- 0
B_C <- PARAMS[1:p_C]
B_x <- PARAMS[(p_C + 1):(p_C + numcol)]
IC <- PARAMS[(p_C + numcol + 1)]
sig <- PARAMS[(p_C + numcol + 2)]
df <- PARAMS[(p_C + numcol + 3)]
}
else if (p_C == 0 && p_NC == 0){
B_NC <- 0
B_C <- 0
B_x <- PARAMS[(p_C + 3):(p_C + 2 + numcol)]
IC <- PARAMS[(p_C + numcol + 3)]
sig <- PARAMS[(p_C + numcol + 4)]
df <- PARAMS[(p_C + numcol + 5)]
}
}
else{
numcol <- 0
B_x <- 0
if (p_C > 0 && p_NC > 0){
B_C <- PARAMS[1:p_C]
B_NC <- PARAMS[(p_C+1):(p_C + p_NC)]
IC <- PARAMS[(p_C + p_NC + 1)]
sig <- PARAMS[(p_C + p_NC + 2)]
df <- PARAMS[(p_C + p_NC + 3)]
}
else if (p_NC > 0 && p_C == 0){
B_C <- 0
B_NC <- PARAMS[1:p_NC]
IC <- PARAMS[(p_NC + 1)]
sig <- PARAMS[(p_NC + 2)]
df <- PARAMS[(p_NC + 3)]
}
else if (p_C > 0 && p_NC == 0){
B_NC <- 0
B_C <- PARAMS[1:p_C]
IC <- PARAMS[(p_C + 1)]
sig <- PARAMS[(p_C + 2)]
df <- PARAMS[(p_C + 3)]
}
else if (p_C == 0 && p_NC == 0){
B_NC <- 0
B_C <- 0
IC <- PARAMS[1]
sig <- PARAMS[2]
df <- PARAMS[3]
}
}
ZC1 <- y[(p_C+1):length(y)]
ZC1 <- fBasics::vec(ZC1)
ZC2 <- regressor.matrix(y,"not",p_C)
if (p_C == 1){
ZC2 <- fBasics::vec(ZC2)
}
if (p_C > 0){
V <- ZC1 - ZC2 %*% B_C
}
else{
V <- ZC1
}
U <- rev(V)
U <- fBasics::vec(U)
ZNC1 <- U[(p_NC + 1):length(U)]
ZNC1 <- fBasics::vec(ZNC1)
ZNC2 <- regressor.matrix(U,"not",p_NC)
if(numcol > 1){
for (i in 1:numcol){
x[,i] <- rev(x[,i])
}
}
else{
x <- rev(x)
}
if(length(x) > 1){
if (numcol > 1 ){
x <- x[(p_NC +1):length(U),]
}
else{
x <- x[(p_NC +1):length(U)]
x <- fBasics::vec(x)
}
}
else{
x <- "not"
}
if (p_NC == 1){
ZNC2 <- fBasics::vec(ZNC2)
}
if (length(x) > 1){
if (p_NC > 0){
E <- rev(ZNC1 - (ZNC2 %*% B_NC) - IC - (x %*% B_x))
}
else{
E <- rev(ZNC1 - IC - (x %*% B_x))
}
}
else{
if (p_NC > 0){
E <- rev(ZNC1 - (ZNC2 %*% B_NC) - IC)
}
else{
E <- rev(ZNC1 - IC)
}
}
se <- sqrt(diag(solve(optimization_results$hessian)))
se.dist <- se[(length(se)-1):length(se)]
se.dist <- rev(se.dist)
return(list(coef.c = B_C, coef.nc = B_NC, coef.exo = B_x, coef.int = IC, scale = sig,df = df,residuals = E, se.dist = se.dist))
}
#' @title Asymptotic inference for the MARX function
#' @description This function allows you to calculate standard errors and confidence intervals for parameters of the MARX model.
#' @param y Data vector of time series observations.
#' @param x Matrix of data (every column represents one time series). Specify NULL or "not" if not wanted.
#' @param B_C Estimated causal parameters of the MARX.
#' @param B_NC Estimated noncausal parameters of the MARX.
#' @param B_x Estimated parameters of the exogenous variables in the MARX.
#' @param IC Estimated intercept.
#' @param sig Estimated scale parameter of the assumed underlying Student-t distribution of the residuals.
#' @param df Estimated degrees of freedom of the assumed underlying Student-t distribution of the residuals.
#' @param sig_level Significance level for the construction of inference.
#' @keywords inference
#' @keywords causal-noncausal
#' @return \item{CI.c}{Confidence intervals for causal parameters.}
#' @return \item{CI.nc}{Confidence intervals for noncausal parameters.}
#' @return \item{CI.exo}{Confidence intervals for exogenous parameters.}
#' @return \item{CI.int}{Confidence interval for intercept.}
#' @return \item{se.c}{Standard errors of causal parameters.}
#' @return \item{se.nc}{Standard errors of noncausal parameters.}
#' @return \item{se.exo}{Standard errors of exogenous parameters.}
#' @return \item{se.int}{Standard error of intercept.}
#' @author Sean Telg
#' @export
#' @examples
#' data <- sim.marx(c('t',1,1), c('t',1,1),100,0.5,0.4,0.3)
#' y <- data$y
#' x <- data$x
#' res <- marx.t(y,x,1,1)
#' inference(y,x,res$coef.c,res$coef.nc,res$coef.exo,res$coef.int,res$scale,res$df,0.05)
inference <- function(y,x,B_C,B_NC,B_x,IC,sig,df,sig_level){
if (is.null(x)){
x <- "not"
}
p_C = length(B_C)
p_NC = length(B_NC)
p_x = length(B_x)
if (length(B_C) == 1 && B_C == 0){
p_C = 0
}
if (length(B_NC) == 1 && B_NC == 0){
p_NC = 0
}
obser = length(y) - p_C - p_NC
regressor_C <- regressor.matrix(y,"not",p_C)
regressand_C <- fBasics::vec(y[(p_C+1):length(y)])
if (p_C == 1){
regressor_C <- fBasics::vec(regressor_C)
}
regressor_NC <- regressor.matrix(fBasics::vec(rev(y)),"not",p_NC)
if (p_NC > 1){
for(i in 1:NCOL(regressor_NC)){
regressor_NC[,i] <- fBasics::vec(rev(regressor_NC[,i]))
}
}
if (p_NC == 1){
regressor_NC <- fBasics::vec(rev(fBasics::vec(regressor_NC)))
}
#if (p_NC > 1){
# for (i in 1:length(regressor_NC[1,])){
# regressor_NC[,i] <- rev(regressor_NC[,i])
# }
#}
regressand_NC <- fBasics::vec(y[1:(length(y)-p_NC)])
## Causal part:
if (p_C > 0){
if (p_NC > 0){
U <- regressand_NC - regressor_NC %*% B_NC
}
else{
U <- regressand_NC
}
U_regressor_C <- regressor.matrix(U,"not",p_C)
Gam_C <- (t(U_regressor_C) %*% U_regressor_C)/obser
Cov_C <- (df+3)/(df+1)*sig^2*solve(Gam_C)
std_c <- matrix(data=NA,nrow=length(Cov_C[,1]), ncol=1)
BC <- matrix(data=NA,nrow=length(Cov_C[,1]), ncol=2)
for(i in 1:length(Cov_C[,1])){
std_c[i] = sqrt(Cov_C[i,i]/obser)
BC[i,] = c((B_C[i] - stats::qnorm(1-(sig_level/2))*std_c[i]),(B_C[i] + stats::qnorm(1-(sig_level/2))*std_c[i]))
}
}
else if (p_C == 0){
BC = 0
std_c = 0
}
## Noncausal part:
if (p_NC > 0){
if (p_C > 0){
V <- regressand_C - regressor_C %*% B_C
}
else{
V <- regressand_C
}
V_regressor_NC <- regressor.matrix(fBasics::vec(rev(V)),"not",p_NC)
if (p_NC > 1){
for(i in 1:NCOL(V_regressor_NC)){
V_regressor_NC[,i] <- fBasics::vec(rev(V_regressor_NC[,i]))
}
}
if (p_NC == 1){
V_regressor_NC <- fBasics::vec(rev(fBasics::vec(V_regressor_NC)))
}
#for (i in 1:length(V_regressor_NC[1,])){
# V_regressor_NC[,i] <- rev(V_regressor_NC[,i])
#}
Gam_NC <- (1/obser)* (t(V_regressor_NC) %*% V_regressor_NC)
Cov_NC <- (df+3)/(df+1)*sig^2*solve(Gam_NC)
std_nc <- matrix(data=NA,nrow=length(Cov_NC[,1]), ncol=1)
BNC <- matrix(data=NA,nrow=length(Cov_NC[,1]), ncol=2)
for(i in 1:length(Cov_NC[,1])){
std_nc[i] = sqrt(Cov_NC[i,i]/obser)
BNC[i,] = c((B_NC[i] - stats::qnorm(1-(sig_level/2))*std_nc[i]),(B_NC[i] + stats::qnorm(1-(sig_level/2))*std_nc[i]))
}
}
else if (p_NC == 0){
BNC = 0
std_nc = 0
}
## Intercept
Cov_IC <- (df+3)/(df+1)*sig^2
std_ic <- sqrt(Cov_IC/obser)
BIC <- c((IC - stats::qnorm(1-(sig_level/2))*std_ic),(IC + stats::qnorm(1-(sig_level/2))*std_ic))
## Coefficients of exogenous variables
if (length(x) > 1){
Gam_x <- (1/obser) * (t(x) %*% x)
Cov_x <- (df+3)/(df+1)*sig^2*solve(Gam_x)
std_x <- matrix(data=NA,nrow=length(Cov_x[,1]), ncol=1)
Bx <- matrix(data=NA,nrow=length(Cov_x[,1]), ncol=2)
for(i in 1:length(Cov_x[,1])){
std_x[i] = sqrt(Cov_x[i,i]/obser)
Bx[i,] = c((B_x[i] - stats::qnorm(1-(sig_level/2))*std_x[i]),(B_x[i] + stats::qnorm(1-(sig_level/2))*std_x[i]))
}
}
else{
Bx <- 0
std_x <- 0
}
return(list(CI.c = BC, CI.nc = BNC, CI.exo = Bx, CI.int = BIC, se.c = std_c, se.nc = std_nc, se.exo = std_x, se.int = std_ic))
}
#' @title The value of the t-log-likelihood for MARX function
#' @description This function allows you to determine the value of the t-log-likelihood for the MARX model.
#' @param params List of parameters.
#' @param y Data vector of time series observations.
#' @param x Matrix of data (every column represents one time series). Specify NULL or "not" if not wanted.
#' @param p_C Number of lags.
#' @param p_NC Number of leads.
#' @keywords optimization
#' @return \item{neg.loglikelihood}{Minus the loglikelihood.}
#' @author Sean Telg
#' @export
#' @examples
#' data <- sim.marx(c('t',1,1), c('t',1,1),100,0.5,0.4,0.3)
#' y <- data$y
#' x <- data$x
#' p_C <- 1
#' p_NC <- 1
#' params <- c(0.5,0.4,0.3,0,1,1)
#' ll.max(params,y,x,p_C,p_NC)
ll.max <- function(params,y,x,p_C,p_NC){
if (is.null(x)){
x <- "not"
}
y <- fBasics::vec(y)
if (length(x) > 1){
colnum <- NCOL(x)
if (p_C > 0 && p_NC > 0){
BC1 <- params[1:p_C]
BNC1 <- params[(p_C+1):(p_C + p_NC)]
Bx1 <- params[(p_C+ p_NC + 1):(p_C + p_NC + colnum)]
IC1 <- params[(p_C + p_NC + colnum + 1)]
sig1 <- params[(p_C + p_NC + colnum + 2)]
df1 <- params[(p_C + p_NC + colnum + 3)]
}
else if (p_NC > 0 && p_C == 0){
BC1 <- 0
BNC1 <- params[1:p_NC]
Bx1 <- params[(p_NC+1):(p_NC + colnum)]
IC1 <- params[(p_NC + colnum + 1)]
sig1 <- params[(p_NC + colnum + 2)]
df1 <- params[(p_NC + colnum + 3)]
}
else if (p_C > 0 && p_NC == 0){
BNC1 <- 0
BC1 <- params[1:p_C]
Bx1 <- params[(p_C + 1): (p_C + colnum)]
IC1 <- params[(p_C + colnum + 1)]
sig1 <- params[(p_C + colnum + 2)]
df1 <- params[(p_C + colnum + 3)]
}
else if (p_C == 0 && p_NC == 0){
BNC1 <- 0
BC1 <- 0
Bx1 <- params[(1:colnum)]
IC1 <- params[(colnum + 1)]
sig1 <- params[(colnum + 2)]
df1 <- params[(colnum + 3)]
}
}
else{
colnum <- 0
if (p_C > 0 && p_NC > 0){
BC1 <- params[1:p_C]
BNC1 <- params[(p_C+1):(p_C + p_NC)]
IC1 <- params[(p_C + p_NC + 1)]
sig1 <- params[(p_C + p_NC + 2)]
df1 <- params[(p_C + p_NC + 3)]
}
else if (p_NC > 0 && p_C == 0){
BC1 <- 0
BNC1 <- params[1:p_NC]
IC1 <- params[(p_NC + 1)]
sig1 <- params[(p_NC + 2)]
df1 <- params[(p_NC + 3)]
}
else if (p_C > 0 && p_NC == 0){
BNC1 <- 0
BC1 <- params[1:p_C]
IC1 <- params[(p_C + 1)]
sig1 <- params[(p_C + 2)]
df1 <- params[(p_C + 3)]
}
else if (p_C == 0 && p_NC == 0){
BNC1 <- 0
BC1 <- 0
IC1 <- params[1]
sig1 <- params[2]
df1 <- params[3]
}
}
ZC1 <- y[(p_C+1):length(y)]
ZC1 <- fBasics::vec(ZC1)
ZC2 <- regressor.matrix(y,"not",p_C)
if (p_C == 1){
ZC2 <- fBasics::vec(ZC2)
}
if (p_C > 0){
V <- ZC1 - (ZC2 %*% BC1)
}
else{
V <- ZC1
}
U <- rev(V)
U <- fBasics::vec(U)
ZNC1 <- U[(p_NC + 1):length(U)]
ZNC1 <- fBasics::vec(ZNC1)
ZNC2 <- regressor.matrix(U,"not",p_NC)
if(colnum > 1){
for (i in 1:colnum){
x[,i] <- rev(x[,i])
}
}
else{
x <- rev(x)
}
if (length(x) > 1){
if (colnum > 1){
x <- x[(p_NC +1):length(U),]
}
else{
x <- x[(p_NC + 1):length(U)]
x <- fBasics::vec(x)
}
}
else{
x = "not"
}
if (p_NC == 1){
ZNC2 <- fBasics::vec(ZNC2)
}
if (length(x) > 1){
if (p_NC > 0){
E <- rev(ZNC1 - (ZNC2 %*% BNC1) - IC1 - (x %*% Bx1))
}
else{
E <- rev(ZNC1 - IC1 - (x %*% Bx1))
}
}
else{
if (p_NC > 0){
E <- rev(ZNC1 - (ZNC2 %*% BNC1) - IC1)
}
else{
E <- rev(ZNC1 - IC1)
}
}
n <- length(E)
loglik_eval <- -(n*lgamma((df1+1)/2) - n*log(sqrt(df1*pi*sig1^2)) - n*lgamma(df1/2) - ((df1+1)/2)*log(1+(E/sig1)^2/df1) %*% matlab::ones(n,1))
return(neg.loglikelihood = loglik_eval)
}
#' @title The lag-lead model selection for MARX function
#' @description This function allows you to determine the MARX model (for p = r + s) that maximizes the t-log-likelihood.
#' @param y Data vector of time series observations.
#' @param x Matrix of data (every column represents one time series). Specify NULL or "not" if not wanted.
#' @param p_pseudo Number of autoregressive terms to be included in the pseudo-causal model.
#' @keywords selection
#' @keywords causal-noncausal
#' @return \item{p.C}{The number of lags selected.}
#' @return \item{p.NC}{The number of leads selected.}
#' @return \item{loglikelihood}{The value of the loglikelihood for all models with p = r + s.}
#' @author Sean Telg
#' @export
#' @examples
#' data <- sim.marx(c('t',3,1), c('t',3,1),100,0.5,0.4,0.3)
#' selection.lag.lead(data$y,data$x,2)
selection.lag.lead <- function(y,x,p_pseudo){
if (is.null(x)){
x <- "not"
}
P_C <- seq(length=(p_pseudo+1), from=0, by=1)
P_C <- fBasics::vec(P_C)
P_NC <- rev(P_C)
P_NC <- fBasics::vec(P_NC)
n <- length(y) - p_pseudo
loglik <- c()
for (i in 1:(p_pseudo+1)){
marx.t_results <- marx.t(y,x,P_C[i],P_NC[i]);
sig <- marx.t_results[[5]]
df <- marx.t_results[[6]]
E <- marx.t_results[[7]]
loglik[i] <- (n*lgamma((df+1)/2) - n*log(sqrt(df*pi*sig^2)) - n*lgamma(df/2) - ((df+1)/2)*log(1+(E/sig)^2/df) %*% matlab::ones(n,1))
}
maxloglik <- which.max(loglik)
P = cbind(P_C,P_NC)
P = fBasics::vec(P[maxloglik,])
p_C <- P[1]
p_NC <- P[2]
return(list(p.C = p_C, p.NC = p_NC,loglikelihood = rev(loglik)))
}
#' @title The model selection for pseudo-ARX function
#' @description This function allows you to calculate AIC, BIC, HQ for pseudo-ARX models.
#' @param y Data vector of time series observations.
#' @param x Matrix of data (every column represents one time series). Specify NULL or "not" if not wanted.
#' @param p_max Maximum number of autoregressive terms to be included.
#' @keywords selection
#' @keywords pseudo-causal
#' @return \item{bic}{Vector containing values BIC for p=0 up to p_max.}
#' @return \item{aic}{Vector containing values AIC for p=0 up to p_max.}
#' @return \item{hq}{vector containing values HQ for p=0 up to p_max.}
#' @author Sean Telg
#' @export
#' @examples
#' data <- sim.marx(c('t',1,1), c('t',1,1),100,0.5,0.4,0.3)
#' selection.lag(data$y,data$x,8)
selection.lag <- function(y,x,p_max){
if (is.null(x)){
x <- "not"
}
bic_results <- bic(y,x,p_max)
aic_results <- aic(y,x,p_max)
hq_results <- hq(y,x,p_max)
bic_vec <- bic_results[[2]]
colnames(bic_vec) <- paste('p =', 0:p_max)
aic_vec <- aic_results[[2]]
colnames(aic_vec) <- paste('p =', 0:p_max)
hq_vec <- hq_results[[2]]
colnames(hq_vec) <- paste('p =', 0:p_max)
cat('Order Selection Criteria Pseudo Causal Model:', "\n")
cat(' ', "\n")
cat('BAYESIAN INFORMATION CRITERION', "\n")
print(bic_vec)
cat(' ', "\n")
cat('Minimum value attained at p = ')
cat(which.min(bic_vec) -1)
cat(' ', "\n")
cat(' ', "\n")
cat('AKAIKE INFORMATION CRITERION', "\n")
print(aic_vec)
cat(' ', "\n")
cat('Minimum value attained at p = ')
cat(which.min(aic_vec) - 1)
cat(' ', "\n")
cat(' ', "\n")
cat('HANNAN-QUINN INFORMATION CRITERION', "\n")
print(hq_vec)
cat(' ', "\n")
cat('Minimum value attained at p = ')
cat(which.min(hq_vec) - 1)
cat(' ', "\n")
cat(' ', "\n")
return(list(bic = bic_vec, aic = aic_vec, hq = hq_vec))
}
#' @title The Bayesian/Schwarz information criterion (BIC) function
#' @description This function allows you to calculate the Bayesian/Schwarz information criteria (BIC) for ARX models.
#' @param y Data vector of time series observations.
#' @param x Matrix of data (every column represents one time series). Specify NULL or "not" if not wanted.
#' @param p_max Maximum number of autoregressive terms to be included.
#' @keywords selection
#' @return \item{p}{Lag order chosen by BIC.}
#' @return \item{values}{Vector containing values BIc for p = 0 up to p_max.}
#' @author Sean Telg
#' @export
#' @examples
#' data <- sim.marx(c('t',1,1), c('t',1,1),100,0.5,0.4,0.3)
#' bic(data$y, data$x,8)
bic <- function(y,x,p_max){
if (is.null(x)){
x <- "not"
}
y <- fBasics::vec(y)
# y <- y - mean(y)
# n <- length(y) - p_max
if (length(x) > 1){
numcol <- NCOL(x)
}
else{
numcol = 0
}
crit <- matrix(data=NA, nrow=(p_max+1), ncol=1)
for (p in 0:p_max){
arx.ls_results <- arx.ls(fBasics::vec(y),x,p)
n <- length(arx.ls_results[[5]])
Cov <- arx.ls_results[[6]]
crit[(p+1)] <- -2*Cov/n + ((log(n))/n)*(p+1+numcol)
}
p_bic <- which.min(crit) - 1
crit <- t(crit)
colnames(crit) <- paste('p =', 0:p_max)
return(list(p = p_bic, values= crit))
}
#' @title The Akaike information criterion (AIC) function
#' @description This function allows you to calculate the Akaike information criteria (AIC) for ARX models.
#' @param y Data vector of time series observations.
#' @param x Matrix of data (every column represents one time series). Specify NULL or "not" if not wanted.
#' @param p_max Maximum number of autoregressive terms to be included.
#' @keywords selection
#' @return \item{p}{Lag order chosen by AIC.}
#' @return \item{values}{Vector containing values AIC for p = 0 up to p_max.}
#' @author Sean Telg
#' @export
#' @examples
#' data <- sim.marx(c('t',1,1), c('t',1,1),100,0.5,0.4,0.3)
#' aic(data$y, data$x,8)
aic <- function(y,x,p_max){
if (is.null(x)){
x <- "not"
}
y <- fBasics::vec(y)
# y <- y - mean(y)
#n <- length(y) - p_max
if (length(x) > 1){
numcol <- NCOL(x)
}
else{
numcol = 0
}
crit <- matrix(data=NA, nrow=(p_max+1), ncol=1)
for (p in 0:p_max){
arx.ls_results <- arx.ls(fBasics::vec(y),x,p)
n <- length(arx.ls_results[[5]])
Cov <- arx.ls_results[[6]]
crit[(p+1)] <- -2*Cov/n + (2/n)*(p+1+numcol)
}
p_aic <- which.min(crit) - 1
crit <- t(crit)
colnames(crit) <- paste('p =', 0:p_max)
return(list(p = p_aic,values = crit))
}
#' @title The Hannan-Quinn (HQ) information criterion function
#' @description This function allows you to calculate the Hannan-Quinn (HQ) information criteria for ARX models.
#' @param y Data vector of time series observations.
#' @param x Matrix of data (every column represents one time series). Specify NULL or "not" if not wanted.
#' @param p_max Maximum number of autoregressive terms to be included.
#' @keywords selection
#' @return \item{p}{Lag order chosen by HQ.}
#' @return \item{values}{Vector containing values HQ for p = 0 up to p_max.}
#' @author Sean Telg
#' @export
#' @examples
#' data <- sim.marx(c('t',1,1), c('t',1,1),100,0.5,0.4,0.3)
#' hq(data$y, data$x,8)
hq <- function(y,x,p_max){
if (is.null(x)){
x <- "not"
}
y <- fBasics::vec(y)
# y <- y - mean(y)
# n <- length(y) - p_max
if (length(x) > 1){
numcol <- NCOL(x)
}
else{
numcol = 0
}
crit <- matrix(data=NA, nrow=(p_max+1), ncol=1)
for (p in 0:p_max){
arx.ls_results <- arx.ls(fBasics::vec(y),x,p)
n <- length(arx.ls_results[[5]])
Cov <- arx.ls_results[[6]]
crit[(p+1)] <- -2*Cov/n + ((2*log(log(n)))/n)*(p+1+numcol)
}
p_hq <- which.min(crit) - 1
crit <- t(crit)
colnames(crit) <- paste('p =', 0:p_max)
return(list(p = p_hq, values = crit))
}
#' @title The simulation of MARX processes
#' @description This function allows you to simulate MARX processes based on different underlying distribution.
#' @param dist.eps vector containing the error distribution and its parameters (options: t, normal, stable).
#' @param dist.x vector containing the distribution of x and its parameters (options: t, normal, stable). Specify NULL or "not" if not wanted.
#' @param obs Number of observations for simulated process.
#' @param c_par vector of causal parameters.
#' @param nc_par vector of noncausal parameters.
#' @param exo_par Parameter of the exogenous variable.
#' @keywords simulation
#' @return \item{y}{Simulated data y.}
#' @return \item{x}{Simulated data x (exogenous variable).}
#' @author Sean Telg
#' @export
#' @examples
#' dist.eps <- c('t',1,1) ## t-distributed errors with 1 degree of freedom and scale parameter 1
#' dist.x <- c('normal',0,1) ## standard normally distributed x variable
#' obs <- 100
#' c_par <- c(0.2,0.4)
#' nc_par <- 0.8
#' exo_par <- 0.5
#' sim.marx(dist.eps,dist.x,obs,c_par,nc_par,exo_par) ## Simulates a MARX(2,1,1) process
sim.marx <- function(dist.eps, dist.x, obs, c_par, nc_par, exo_par){
if (is.null(dist.x)){
exo_par = 0
dist.x = 0
}
if (length(dist.x) == 1){
exo_par = 0
}
M <- 5*obs + 300
M.star <- 500 #0.05*M
if (dist.eps[1] == 't'){
eps <- as.numeric(dist.eps[3])*(stats::rt((3*(obs+200)), df = as.numeric(dist.eps[2]), ncp = 0))
} else if (dist.eps[1] == 'cauchy'){
eps <- stats::rcauchy((3*(obs+200)), location= as.numeric(dist.eps[2]), scale=as.numeric(dist.eps[3]))
} else if (dist.eps[1] == 'normal') {
eps <- stats::rnorm((3*(obs+200)), mean = as.numeric(dist.eps[2]), sd=as.numeric(dist.eps[3]))
} else if (dist.eps[1] == 'stable') {
eps <- stabledist::rstable((3*(obs+200)), alpha = as.numeric(dist.eps[2]), beta=as.numeric(dist.eps[3]), gamma=as.numeric(dist.eps[4]), delta=as.numeric(dist.eps[5]))
}
if (dist.x[1] == 't'){
x <- as.numeric(dist.x[3])*(stats::rt((3*(M+200)), df = as.numeric(dist.x[2]), ncp = 0))
} else if (dist.x[1] == 'cauchy'){
x <- stats::rcauchy((3*(M+200)), location= as.numeric(dist.x[2]), scale=as.numeric(dist.x[3]))
} else if (dist.x[1] == 'normal') {
x <- stats::rnorm((3*(M+200)), mean = as.numeric(dist.x[2]), sd=as.numeric(dist.x[3]))
} else if (dist.x[1] == 'stable') {
x <- stabledist::rstable((3*(M+200)), alpha = as.numeric(dist.x[2]), beta=as.numeric(dist.x[3]), gamma=as.numeric(dist.x[4]), delta=as.numeric(dist.x[5]))
}
eps <- eps[(2*(obs+200)):(3*(obs+200))]
r <- length(c_par)
s <- length(nc_par)
if (length(dist.x) != 1){
BC <- c(c_par, rep(0,M.star - length(c_par)))
Phi.matrix <- matrix(data=0, nrow=M.star, ncol=M.star)
Psi <- c()
Psi[1] <- 1
for (i in 1:M.star){
for (j in 1:i){
Phi.matrix[i,j] = Psi[i-j+1]*BC[j]
}
Psi[(i+1)] <- sum(Phi.matrix[i,])
}
## Generating linear filter over x_t "similar to MA representation for AR"
MA_rep <- matrix(data=NA, nrow=(obs+200), ncol=M.star)
for (i in 1:(obs+200)){
for (j in 1:M.star){
MA_rep[i,j] <- Psi[j]*x[M-(j-1)-(i-1)]
}
}
exo_var <- c()
exo_var <- MA_rep[,1]
MA_rep <- rowSums(MA_rep)
#MA_rep <- rev(MA_rep)
}
else{
MA_rep <- matrix(data=0, nrow=(obs+200),ncol=M.star)
exo_var <- c()
}
u <- matrix(data=NA, nrow=(obs+100+s), ncol=1)
U.mat <- u[(201+obs):(200+obs+s)] <- rep(1,s)
for (i in (obs+200):1){
if (length(nc_par) == 1){
AR <- nc_par * U.mat
}
else{
AR <- sum(t(nc_par*U.mat))
}
if (i == (obs+200)){AR <- 0}
u[i] <- AR + eps[i]
if (s == 1){
U.mat <- u[i]
}
else{
U.mat <- c(u[i], t(U.mat[1:(length(U.mat)-1)]))
U.mat <- t(U.mat)
}
}
y <- matrix(data=0, nrow=(obs+r+200), ncol=1)
Y <- y[1:r] <- rep(0,r)
for (i in 1:(obs+200)){
if (r == 1){
AR2 <- c_par * Y
}
else{
AR2 <- sum(t(c_par*Y))
}
if (i == 1){AR2 <- 0}
y[i+r] <- AR2 + u[i] + exo_par * MA_rep[i+r]
if (r == 1){
Y <- y[i+r]
}
else{
Y <- c(y[i+r], t(Y[1:(length(Y)-1)]))
Y <- t(Y)
}
}
u <- u[1:obs]
y <- y[(101+r):(obs+100+r)]
exo_var <- exo_var[(101+r):(obs+100+r)]
return(list(y = y, x = exo_var))
}
#' @title Companion form function
#' @description This function allows you to compute a companion form matrix in order to check the stability of causal and noncausal part of the ARX model.
#' @param pol Coefficient vector. If polynomial is 1 - ax - bx^2, coefficient vector is c(a, b).
#' @keywords stability, stationarity
#' @return \item{C}{Companion matrix C.}
#' @author Sean Telg
#' @export
#' @examples
#' pol <- c(0.3,0.4)
#' C <- companion.form(pol)
companion.form <- function(pol){
r <- length(pol)
b <- cbind(diag((r-1)),rep(0,(r-1)))
C <- rbind(t(pol),b)
return(C)
}
#' @title Coefficients of the moving average representation function
#' @description This function allows you to invert a polynomial (either the causal or the noncausal one) and output the corresponding coefficients of the moving average representation.
#' @param pol Coefficient vector. If polynomial is 1 - ax - bx^2, coefficient vector is c(a, b).
#' @param M Truncation value M (how many MA coefficients should be computed?).
#' @keywords stability, stationarity
#' @return \item{psi}{Vector containing coefficients of the moving average representation.}
#' @author Sean Telg
#' @export
#' @examples
#' pol <- c(0.3,0.4)
#' psi <- companion.form(pol)
compute.MA <- function(pol,M){
r <- length(pol)
str <- c(rep(0,(r-1)),1)
psi <- c(1)
for (j in 1:M){
psi[j+1] <- t(pol) %*% rev(str)
if (r > 1){
str <- c(str[2:length(str)],psi[j+1])
}
else{
str <- psi[j+1]
}
}
return(psi)
}
#' @title Forecasting function for the MARX model
#' @description This function allows you to forecast with the mixed causal-noncausal model with possibly exogenous regressors.
#' @param y Data vector y.
#' @param X (optional) Matrix with data (column represent a series).
#' @param p_C Number of lags (causal order).
#' @param p_NC Number of leads (noncausal order).
#' @param X.for (optional) Matrix with forecasted values for X (column represents series).
#' @param h Forecast horizon h.
#' @param M (optional) Truncation value M for MA representation. Default value: 50.
#' @param N (optional) Number of simulations to forecast noncausal component. Default: 10,000.
#' @keywords forecasting
#' @return \item{y.for}{Vector containing forecasted values for y.}
#' @author Sean Telg
#' @export
#' @examples
#' ## Forecasting MAR(0,1) model 4-periods ahead for lnbev (from dataset)
#' data <- MARX::dataset[,2]
#' y.for <- forecast.marx(y=data, p_C=0, p_NC=1, h=4, M=50, N=1000)
forecast.marx <- function(y,X,p_C,p_NC,X.for,h,M,N){
set.seed(9999)
if (missing(X) == TRUE){
X = NULL
}
if (missing(N) == TRUE){
N = 10000
}
object <- mixed(y,X,p_C,p_NC)
obs <- length(y)
## Check whether there are exogenous variables and whether truncation M is known
if (missing(X.for) == TRUE && missing(M) == TRUE){
X.for = NULL
M = 50
}
else if(missing(X.for) == TRUE && missing(M) == FALSE){
X.for = NULL
M = M
}
else if(missing(X.for) == FALSE && missing(M) == TRUE){
if (NCOL(X.for) == 1){
if(is.null(X.for) == TRUE){
M = 50
}
else{
M = length(X.for)
}
}
else{
M = length(X.for[,1])
}
}
else if(missing(X.for) == FALSE && missing(M) == FALSE){
if (NCOL(X.for) == 1){
if(is.null(X.for) == TRUE){
M = M
}
else{
M = min(length(X.for), M)
}
}
else{
M = min(length(X.for[,1]),M)
}
}
coef.caus <- c()
if (object$order[1] == 0){
r = 1
coef.caus <- object$coefficients[(r+1)]
}
else{
r = object$order[1]
coef.caus <- object$coefficients[2:(r+1)]
}
coef.noncaus <- c()
if (object$order[2] == 0){
s = 1
coef.noncaus <- object$coefficients[(r+1+s)]
}
else{
s = object$order[2]
coef.noncaus <- object$coefficients[(r+2):(r+1+s)]
}
coef.exo <- c()
if (object$order[3] == 0){
q = 1
coef.exo <- object$coefficients[(r+1+s+q)]
}
else{
q = object$order[3]
coef.exo <- object$coefficients[(r+1+s+1):(r+s+1+q)]
}
## Simulate future epsilon and use forecasted X
hve <- c()
hve2 <- matrix(data=0, nrow=N,ncol=h)
for (iter in 1:N){
eps.sim <- object$coefficients["scale",]*stats::rt(M,object$coefficients["df",])
z2 <- c()
for (i in 1:M){
if(is.null(X.for) == TRUE){
z2[i] <- eps.sim[i]
}
else{
if(NCOL(X.for) > 1){
z2[i] <- eps.sim[i] + coef.exo %*% t(X.for[i,])
}
else{
z2[i] <- eps.sim[i] + coef.exo * X.for[i]
}
}
}
## Compute filtered values u = phi(L)y and moving average values
phi <- c(1,coef.caus)
u <- c()
for (i in (r+1):obs){
u[i] <- phi %*% y[i:(i-r)]
}
w <- c(u[(obs-s+1):obs],z2)
C <- matrix(data=0, nrow=(M+s), ncol=(M+s))
C[1,] <- compute.MA(coef.noncaus,(M+s-1))
if (s > 1){
for (i in 2:s){
C[i,] <- c(0, C[(i-1),1:(length(C[(i-1),])-1)])
}
}
for (i in (s+1):(M+s)){
C[i,] <- c(rep(0,(i-1)),1,rep(0,(M+s-i)))
}
D = solve(C)
e <- D %*% w
h1 <- c()
for (i in 1:s){
h1[i] <- metRology::dt.scaled(e[i], df=object$coefficients["df",], sd=object$coefficients["scale",])
}
hve[iter] = prod(h1)
for (j in 1:h){
mov.av <- C[1,1:(M-j+1)] %*% z2[j:M]
hve2[iter,j] <- mov.av * hve[iter]
}
}
y.star <- y[(obs-r+1):obs]
y.for <- c()
exp <- c()
for (j in 1:h){
exp[j] = ((1/N)*sum(hve2[,j]))/((1/N)*sum(hve))
if(length(coef.caus) == 1){
y.for[j] <- object$coefficients[1]/(1-sum(coef.noncaus)) + coef.caus * y.star + exp[j]
}
else{
y.for[j] <- object$coefficients[1]/(1-sum(coef.noncaus)) + t(coef.caus) %*% y.star + exp[j]
}
y.star <- c(y.for[j], y.star[1:(length(y.star)-1)])
}
return(y.for)
}
mixed.combine <- function(y,x,p_C,p_NC){
q <- NCOL(x)
if (is.null(x) == TRUE){
x <- "not"
q <- 0
}
sig_level = 0.05
est <- marx.t(y,x,p_C,p_NC)
residuals <- est$residuals
fitted.values <- y[(1+p_C):(length(y)-p_NC)] - residuals
inf <- inference(y,x,est$coef.c,est$coef.nc,est$coef.exo,est$coef.int,est$scale,est$df,sig_level)
coefficients <- c(est$coef.int,est$coef.c,est$coef.nc,est$coef.exo, est$df, est$scale)
coefficients <- fBasics::vec(coefficients)
if(length(x) > 1){
if (p_C > 0 && p_NC > 0){
rownames(coefficients) <- c('int', paste('lag', 1:p_C), paste('lead', 1:p_NC), paste('exo', 1:NCOL(x)), 'df', 'scale')
}
else if (p_C > 0 && p_NC == 0){
rownames(coefficients) <- c('int', paste('lag', 1:p_C), 'lead', paste('exo', 1:NCOL(x)), 'df', 'scale')
}
else if (p_C == 0 && p_NC > 0){
rownames(coefficients) <- c('int', 'lag', paste('lead', 1:p_NC), paste('exo', 1:NCOL(x)), 'df', 'scale')
}
else if (p_C == 0 && p_NC == 0){
rownames(coefficients) <- c('int', 'lag', 'lead', paste('exo', 1:NCOL(x)), 'df', 'scale')
}
}
else{
if (p_C > 0 && p_NC > 0){
rownames(coefficients) <- c('int', paste('lag', 1:p_C), paste('lead', 1:p_NC), 'exo', 'df', 'scale')
}
else if (p_C > 0 && p_NC == 0){
rownames(coefficients) <- c('int', paste('lag', 1:p_C), 'lead', 'exo', 'df', 'scale')
}
else if (p_C == 0 && p_NC > 0){
rownames(coefficients) <- c('int', 'lag', paste('lead', 1:p_NC), 'exo', 'df', 'scale')
}
else if (p_C == 0 && p_NC == 0){
rownames(coefficients) <- c('int', 'lag', 'lead', 'exo', 'df', 'scale')
}
}
se <- c(inf$se.int, inf$se.c, inf$se.nc, inf$se.exo, est$se.dist[1], est$se.dist[2])
if (length(x) > 1){
degree <- NROW(est$residuals) - (p_C + p_NC + NCOL(x) + 1)
}
else{
degree <- NROW(est$residuals) - (p_C + p_NC + 1)
}
return(list(coefficients = coefficients, se = se, df.residual = degree, residuals= residuals, fitted.values = fitted.values, order=c(p_C,p_NC,q)))
}
#' @title The pseudo-causal model function
#' @description This function allows you to estimate pseudo-causal ARX models by OLS (compatible with most functions in lm() class).
#' @aliases pseudo
#' pseudo.default
#' print.pseudo
#' summary.pseudo
#' @param y Data vector of time series observations.
#' @param x Matrix of data (every column represents one time series). Specify NULL or "not" if not wanted.
#' @param p Number of lags to be included.
#' @param object An object of the class "pseudo"
#' @param ... Other arguments
#' @keywords pseudo-causal
#' @keywords estimation
#' @examples
#' data <- sim.marx(c('t',1,1), c('t',1,1),100,0.5,0.4,0.3)
#' object <- pseudo(data$y, data$x, 2)
#' class(object) <- "pseudo"
#' summary(object)
#' @rdname pseudo
#' @export pseudo
pseudo <- function(y,x,p){ UseMethod("pseudo")}
#' @return An object of class \code{"pseudo"} is a list containing the following components:
#' @return \item{coefficients}{Vector of estimated coefficients.}
#' @return \item{coef.auto}{Vector of estimated autoregressive parameters.}
#' @return \item{coef.exo}{Vector of estimated exogenous parameters.}
#' @return \item{mse}{Mean squared error.}
#' @return \item{residuals}{Residuals.}
#' @return \item{loglikelihood}{Value of the loglikelihood.}
#' @return \item{fitted.values}{Fitted values.}
#' @return \item{df}{Degrees of freedom.}
#' @return \item{vcov}{Variance-covariance matrix of residuals.}
#'
#' @rdname pseudo
#' @method pseudo default
#' @S3method pseudo default
pseudo.default <- function(y,x,p){
if (is.null(x)){
x <- "not"
}
est <- arx.ls(y,x,p)
est$call <- match.call();
class(est) <- "pseudo"
est
}
#' @rdname pseudo
#' @method print pseudo
#' @S3method print pseudo
print.pseudo <- function(x,...){
cat("Call:\n")
print(x$call)
cat("\n Coefficients:\n")
print(x$coefficients)
}
#' @rdname pseudo
#' @method summary pseudo
#' @S3method summary pseudo
summary.pseudo <- function(object,...){
se <- sqrt(diag(object$vcov))
tval <- stats::coef(object) / se
TAB <- cbind(Estimate = stats::coef(object),
StdErr = se,
t.value = tval,
p.value = 2*stats::pt(-abs(tval),df=object$df))
colnames(TAB) <- c("Estimate", "Std.Err", "t value", "p value")
res <- list(call=object$call, coefficients = TAB)
class(res) <- "summary.pseudo"
res
}
pseudo <- function(y,x,p){ UseMethod("pseudo")}
pseudo.default <- function(y,x,p){
if (is.null(x)){
x <- "not"
}
est <- arx.ls(y,x,p)
est$call <- match.call();
class(est) <- "pseudo"
est
}
print.pseudo <- function(x,...){
cat("Call:\n")
print(x$call)
cat("\n Coefficients:\n")
print(x$coefficients)
}
summary.pseudo <- function(object,...){
se <- sqrt(diag(object$vcov))
tval <- stats::coef(object) / se
TAB <- cbind(Estimate = stats::coef(object),
StdErr = se,
t.value = tval,
p.value = 2*stats::pt(-abs(tval),df=object$df))
colnames(TAB) <- c("Estimate", "Std.Err", "t value", "p value")
res <- list(call=object$call, coefficients = TAB)
class(res) <- "summary.pseudo"
res
}
#' @title The MARX estimation function
#' @description This function allows you to estimate mixed causal-noncausal MARX models by t-MLE (compatible with most functions in lm() class).
#' @aliases mixed
#' mixed.default
#' print.mixed
#' summary.mixed
#' @param y Data vector of time series observations.
#' @param x Matrix of data (every column represents one time series). Specify NULL or "not" if not wanted.
#' @param p_C Number of lags to be included.
#' @param p_NC Number of leads to be included.
#' @param object An object of the class "mixed".
#' @param ... Other parameters.
#' @keywords causal-noncausal
#' @keywords estimation
#' @examples
#' data <- sim.marx(c('t',1,1), c('t',1,1),100,0.5,0.4,0.3)
#' object <- mixed(data$y, data$x, 1, 1)
#' class(object) <- "mixed"
#' summary(object)
#' @rdname mixed
#' @export mixed
mixed <- function(y,x,p_C,p_NC){ UseMethod("mixed")}
#' @return An object of class \code{"mixed"} is a list containing the following components:
#' @return \item{coefficients}{Vector of estimated coefficients.}
#' @return \item{se}{Standard errors of estimated coefficients.}
#' @return \item{df.residual}{Degrees of freedom residuals.}
#' @return \item{residuals}{Residuals.}
#' @return \item{fitted.values}{Fitted values.}
#' @return \item{order}{Vector containing (r,s,q), i.e. causal order r, noncausal order s, number of exogenous regressors q.}
#' @rdname mixed
#' @method mixed default
#' @S3method mixed default
mixed.default <- function(y,x,p_C,p_NC){
est <- mixed.combine(y,x,p_C,p_NC);
est$call <- match.call();
class(est) <- "mixed"
est
}
#' @rdname mixed
#' @method print mixed
#' @s3method print mixed
print.mixed <- function(x,...){
cat("Call:\n")
print(x$call)
cat("\n Coefficients:\n")
print(x$coefficients)
}
#' @rdname mixed
#' @method summary mixed
#' @S3method summary mixed
summary.mixed <- function(object,...){
tval <- stats::coef(object) / object$se
TAB <- cbind(Estimate = stats::coef(object),
StdErr = object$se,
t.value = tval,
p.value = 2*stats::pt(-abs(tval),df=object$df.residual))
colnames(TAB) <- c("Estimate", "Std.Err", "t value", "p value")
res <- list(call=object$call, coefficients = TAB)
class(res) <- "summary.mixed"
res
}
mixed <- function(y,x,p_C,p_NC){ UseMethod("mixed")}
mixed.default <- function(y,x,p_C,p_NC){
est <- mixed.combine(y,x,p_C,p_NC);
est$call <- match.call();
class(est) <- "mixed"
est
}
print.mixed <- function(x,...){
cat("Call:\n")
print(x$call)
cat("\n Coefficients:\n")
print(x$coefficients)
}
summary.mixed <- function(object,...){
tval <- stats::coef(object) / object$se
TAB <- cbind(Estimate = stats::coef(object),
StdErr = object$se,
t.value = tval,
p.value = 2*stats::pt(-abs(tval),df=object$df.residual))
colnames(TAB) <- c("Estimate", "Std.Err", "t value", "p value")
res <- list(call=object$call, coefficients = TAB)
class(res) <- "summary.mixed"
res
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.