R/makeMatrices.R

In pli2016/forecast: Forecasting Functions for Time Series and Linear Models

# These functions make the w, F, x and g matrices
# 
# 
# Author: srazbash
###############################################################################

makeTBATSFMatrix <- function(alpha, beta=NULL, small.phi=NULL, seasonal.periods=NULL, k.vector=NULL, gamma.bold.matrix=NULL, ar.coefs=NULL, ma.coefs=NULL) {
	
	#1. Alpha Row
	F <- matrix(1,nrow=1,ncol=1)
	if(!is.null(beta)) {
		F <- cbind(F, matrix(small.phi,nrow=1,ncol=1))
	}
	if(!is.null(seasonal.periods)) {
		tau <- sum(k.vector)*2
		zero.tau <- matrix(0,nrow=1, ncol=tau)
		F <- cbind(F,zero.tau)
	}
	if(!is.null(ar.coefs)) {
		p <- length(ar.coefs)
		ar.coefs <- matrix(ar.coefs, nrow=1, ncol=p)
		alpha.phi <- alpha*ar.coefs
		F <- cbind(F, alpha.phi)
	}
	if(!is.null(ma.coefs)) {
		q <- length(ma.coefs) 
		ma.coefs <- matrix(ma.coefs,nrow=1,ncol=q)
		alpha.theta <- alpha*ma.coefs
		F <- cbind(F,alpha.theta)
	}
	
	#2. Beta Row
	if(!is.null(beta)) {
		beta.row <- matrix(c(0,small.phi),nrow=1,ncol=2)
		if(!is.null(seasonal.periods)) {
			beta.row <- cbind(beta.row, zero.tau)
		}
		if(!is.null(ar.coefs)) {
			beta.phi <- beta*ar.coefs
			beta.row <- cbind(beta.row, beta.phi)
		}
		if(!is.null(ma.coefs)) {
			beta.theta <- beta*ma.coefs
			beta.row <- cbind(beta.row, beta.theta)
		}
		F <- rbind(F, beta.row)
	}
	
	#3. Seasonal Row
	if(!is.null(seasonal.periods)) {
		seasonal.row <- t(zero.tau)
		if(!is.null(beta)) {
			seasonal.row <- cbind(seasonal.row,seasonal.row)
		}
		
		
		#Make the A matrix
		A <- matrix(0,tau,tau)
		last.pos <- 0
		for(i in 1:length(k.vector)) {
			if(seasonal.periods[i] != 2) {
				C <- .Call("makeCIMatrix", k_s = as.integer(k.vector[i]), m_s = as.double(seasonal.periods[i]), PACKAGE = "forecast")
			} else {
				C <- matrix(0,1,1)	
			}
			S <- .Call("makeSIMatrix", k_s = as.integer(k.vector[i]), m_s = as.double(seasonal.periods[i]), PACKAGE = "forecast")
			
			#C <- matrix(0,k.vector[i],k.vector[i])
			#for(j in 1:k.vector[i]) {
			#	l <- round((2*pi*j/seasonal.periods[i]), digits=15)
			#	C[j,j] <- cos(l)
			#}
			#S <- matrix(0,k.vector[i],k.vector[i])
			#for(j in 1:k.vector[i]) {
			#	S[j,j] <- sin(2*pi*j/seasonal.periods[i])
			#}
			#print(C)
			#print(S)
			Ai <- .Call("makeAIMatrix", C_s = C, S_s = S, k_s = as.integer(k.vector[i]), PACKAGE = "forecast")
			A[(last.pos+1):(last.pos+(2*k.vector[i])), (last.pos+1):(last.pos+(2*k.vector[i]))] <- Ai
			last.pos <- last.pos+(2*k.vector[i])
			
				
			
		}
		seasonal.row <- cbind(seasonal.row,A)
		
		if(!is.null(ar.coefs)) {
			B <- t(gamma.bold.matrix) %*% ar.coefs
			seasonal.row <- cbind(seasonal.row, B)
		}
		if(!is.null(ma.coefs)) {
			C <- t(gamma.bold.matrix) %*% ma.coefs
			seasonal.row <- cbind(seasonal.row, C)
		}
		F <- rbind(F, seasonal.row)
	}
	
	#4. AR() Rows
	if(!is.null(ar.coefs)) {
		#p <- length(ar.coefs)
		ar.rows <- matrix(0,nrow=p,ncol=1)
		if(!is.null(beta)) {
			ar.rows <- cbind(ar.rows, ar.rows)
		}
		if(!is.null(seasonal.periods)) {
			ar.seasonal.zeros <- matrix(0,nrow=p,ncol=tau)
			ar.rows <- cbind(ar.rows, ar.seasonal.zeros)
		}
		ident <- diag((p-1))
		ident <- cbind(ident, matrix(0,nrow=(p-1),ncol=1))
		ar.part <- rbind(ar.coefs,ident)
		ar.rows <- cbind(ar.rows, ar.part)
		
		if(!is.null(ma.coefs)) {
			ma.in.ar <- matrix(0,nrow=p,ncol=q)
			ma.in.ar[1,] <- ma.coefs
			ar.rows <- cbind(ar.rows,ma.in.ar)
		}
		
		F <- rbind(F,ar.rows)
	}
	
	#5. MA() Rows
	if(!is.null(ma.coefs)) {
		ma.rows <- matrix(0,nrow=q,ncol=1)
		if(!is.null(beta)) {
			ma.rows <- cbind(ma.rows,ma.rows)
		}
		if(!is.null(seasonal.periods)) {
			ma.seasonal <- matrix(0, nrow=q, ncol=tau)
			ma.rows <- cbind(ma.rows,ma.seasonal)
		}
		if(!is.null(ar.coefs)) {
			ar.in.ma <- matrix(0, nrow=q, ncol=p)
			ma.rows <- cbind(ma.rows,ar.in.ma)
		}
		ident <- diag((q-1))
		ident <- cbind(ident,matrix(0,nrow=(q-1),ncol=1))
		ma.part <- rbind(matrix(0,nrow=1,ncol=q),ident)
		ma.rows <- cbind(ma.rows,ma.part)
		F <- rbind(F,ma.rows)
	}
	return(F)
}


#makeWMatrix <- function(small.phi=NULL, seasonal.periods=NULL, ar.coefs=NULL, ma.coefs=NULL) {
#
#	the.list <- .Call("makeBATSWMatrix", smallPhi_s = small.phi, sPeriods_s = as.integer(seasonal.periods), arCoefs_s = ar.coefs, maCoefs_s = ma.coefs, PACKAGE = "forecast")
#
#	
#	return(the.list)
#	
#}

#makeGMatrix <- function(alpha, beta=NULL, gamma.vector=NULL, seasonal.periods=NULL, p=0, q=0) {
#	li <- .Call("makeBATSGMatrix", alpha, beta, gamma.vector, as.integer(seasonal.periods), as.integer(p), as.integer(q), PACKAGE="forecast")
#	
#	return(li)
#}


makeFMatrix <- function(alpha, beta=NULL, small.phi=NULL, seasonal.periods=NULL, gamma.bold.matrix=NULL, ar.coefs=NULL, ma.coefs=NULL) {
	
	#1. Alpha Row
	F <- matrix(1,nrow=1,ncol=1)
	if(!is.null(beta)) {
		F <- cbind(F, matrix(small.phi,nrow=1,ncol=1))
	}
	if(!is.null(seasonal.periods)) {
		tau <- sum(seasonal.periods)
		zero.tau <- matrix(0,nrow=1, ncol=tau)
		F <- cbind(F,zero.tau)
	}
	if(!is.null(ar.coefs)) {
		p <- length(ar.coefs)
		ar.coefs <- matrix(ar.coefs, nrow=1, ncol=p)
		alpha.phi <- alpha*ar.coefs
		F <- cbind(F, alpha.phi)
	}
	if(!is.null(ma.coefs)) {
		q <- length(ma.coefs) 
		ma.coefs <- matrix(ma.coefs,nrow=1,ncol=q)
		alpha.theta <- alpha*ma.coefs
		F <- cbind(F,alpha.theta)
	}

	#2. Beta Row
	if(!is.null(beta)) {
		beta.row <- matrix(c(0,small.phi),nrow=1,ncol=2)
		if(!is.null(seasonal.periods)) {
			beta.row <- cbind(beta.row, zero.tau)
		}
		if(!is.null(ar.coefs)) {
			beta.phi <- beta*ar.coefs
			beta.row <- cbind(beta.row, beta.phi)
		}
		if(!is.null(ma.coefs)) {
			beta.theta <- beta*ma.coefs
			beta.row <- cbind(beta.row, beta.theta)
		}
		F <- rbind(F, beta.row)
	}
	
	#3. Seasonal Row
	if(!is.null(seasonal.periods)) {
		seasonal.row <- t(zero.tau)
		if(!is.null(beta)) {
			seasonal.row <- cbind(seasonal.row,seasonal.row)
		}
		
		
		#Make the A matrix
		for(i in seasonal.periods) {
			if(i == seasonal.periods[1]) {
				a.row.one <- matrix(0,nrow=1,ncol=i)
				a.row.one[i] <- 1
				a.row.two <- cbind(diag((i-1)), matrix(0,nrow=(i-1),ncol=1))
				A <- rbind(a.row.one,a.row.two)
			} else {
				old.A.rows <- dim(A)[1]
				old.A.columns <- dim(A)[2]
				a.row.one <- matrix(0,nrow=1,ncol=i)
				a.row.one[i] <- 1
				a.row.two <- cbind(diag((i-1)), matrix(0,nrow=(i-1),ncol=1))
				Ai <- rbind(a.row.one,a.row.two)
				A <- rbind(A, matrix(0, nrow=dim(Ai)[1], ncol=old.A.columns))
				A <- cbind(A, matrix(0,nrow=dim(A)[1],ncol=dim(Ai)[2]))
				A[((old.A.rows+1):(old.A.rows+dim(Ai)[1])),((old.A.columns+1):(old.A.columns+dim(Ai)[2]))] <- Ai
			}
		}
		seasonal.row <- cbind(seasonal.row,A)
		
		if(!is.null(ar.coefs)) {
			B <- t(gamma.bold.matrix) %*% ar.coefs
			seasonal.row <- cbind(seasonal.row, B)
		}
		if(!is.null(ma.coefs)) {
			C <- t(gamma.bold.matrix) %*% ma.coefs
			seasonal.row <- cbind(seasonal.row, C)
		}
		F <- rbind(F, seasonal.row)
	}
	
	#4. AR() Rows
	if(!is.null(ar.coefs)) {
		#p <- length(ar.coefs)
		ar.rows <- matrix(0,nrow=p,ncol=1)
		if(!is.null(beta)) {
			ar.rows <- cbind(ar.rows, ar.rows)
		}
		if(!is.null(seasonal.periods)) {
			ar.seasonal.zeros <- matrix(0,nrow=p,ncol=tau)
			ar.rows <- cbind(ar.rows, ar.seasonal.zeros)
		}
		ident <- diag((p-1))
		ident <- cbind(ident, matrix(0,nrow=(p-1),ncol=1))
		ar.part <- rbind(ar.coefs,ident)
		ar.rows <- cbind(ar.rows, ar.part)
		
		if(!is.null(ma.coefs)) {
			ma.in.ar <- matrix(0,nrow=p,ncol=q)
			ma.in.ar[1,] <- ma.coefs
			ar.rows <- cbind(ar.rows,ma.in.ar)
		}
		
		F <- rbind(F,ar.rows)
	}

	#5. MA() Rows
	if(!is.null(ma.coefs)) {
		ma.rows <- matrix(0,nrow=q,ncol=1)
		if(!is.null(beta)) {
			ma.rows <- cbind(ma.rows,ma.rows)
		}
		if(!is.null(seasonal.periods)) {
			ma.seasonal <- matrix(0, nrow=q, ncol=tau)
			ma.rows <- cbind(ma.rows,ma.seasonal)
		}
		if(!is.null(ar.coefs)) {
			ar.in.ma <- matrix(0, nrow=q, ncol=p)
			ma.rows <- cbind(ma.rows,ar.in.ma)
		}
		ident <- diag((q-1))
		ident <- cbind(ident,matrix(0,nrow=(q-1),ncol=1))
		ma.part <- rbind(matrix(0,nrow=1,ncol=q),ident)
		ma.rows <- cbind(ma.rows,ma.part)
		F <- rbind(F,ma.rows)
	}
	return(F)
}


makeXMatrix <- function(l, b=NULL, s.vector=NULL, d.vector=NULL, epsilon.vector=NULL) {
	x.transpose <- matrix(l, nrow=1, ncol=1)
	if(!is.null(b)) {
		x.transpose <- cbind(x.transpose, matrix(b, nrow=1, ncol=1))
	}
	if(!is.null(s.vector)) {
		x.transpose <- cbind(x.transpose,matrix(s.vector,nrow=1,ncol=length(s.vector)))
	}
	
	if(!is.null(d.vector)) {
		x.transpose <- cbind(x.transpose,matrix(d.vector,nrow=1,ncol=length(d.vector)))
	}
		
	if(!is.null(epsilon.vector)) {
		x.transpose <- cbind(x.transpose,matrix(epsilon.vector,nrow=1,ncol=length(epsilon.vector)))
	}
	
	x <- t(x.transpose)
	return(list(x=x, x.transpose=x.transpose))
}