# R/MS.R In JieGroup/tsms: Modeling Univariate Time Series with Multiple Seasonality

#### Documented in MS

```MS = function(X,W=NULL,U,p_range,q_range,r=1,S=NULL,
blur.out=c(2,2),sar=T,sma=T,
sfourier=F,order=1,
sfactor=F,pred_t=0,
crit='BC',level=0.05,multicore=T){

# step 1: prepare
X = as.matrix(X)

# step 2: suggest seasonality
if (is.null(S)){
S = suggest(X,level)
}
S = unique(S[(S>(max(p_range,q_range)+blur.out[2]))&(S<(U-blur.out[1]))])

# step 3: prune seasonality
S.nat = sort(S)
S.refer = list()
j=1
for (k in 1:length(S.nat)){
if (k==1){
S.refer[[j]] = S.nat[k]
j=j+1
}else{
S.temp = c()
for (i in 1:length(S.refer[[j-1]])){
S.temp = c(S.temp,c((S.refer[[j-1]][i]-sum(blur.out)):
S.refer[[j-1]][i]+sum(blur.out)))
}
if (S.nat[k]%in%S.temp){
S.refer[[j-1]] = c(S.refer[[j-1]],S.nat[k])
}else{
S.refer[[j]] = S.nat[k]
j = j+1
}
}
}
S.pick = list()
for (k in 1:length(S.refer)){
S.temp = c()
for (j in 1:length(S.refer[[k]])){
S.temp = sort(unique(c(S.temp,c((S.refer[[k]][j]-blur.out[1]):
(S.refer[[k]][j]+blur.out[2])))))
}
S.pick[[k]] = S.temp
}
r.nomial = min(r,length(S.pick))

# step 4: merge W.list
W.list = list(exo = W)
if (sfourier|sfactor){
for (i in 1:length(S.pick)){
S = mean(S.pick[[i]])
sea = matrix(nrow=as.numeric(sfourier)*order*2+as.numeric(sfactor)*S,ncol=nrow(X)+pred_t)
if (sfourier){
for (j in 1:order){
sea[j,] = sin(2*pi*c(1:(nrow(X)+pred_t))/(S/j))
sea[j+order,] = cos(2*pi*c(1:(nrow(X)+pred_t))/(S/j))
}
}
if (sfactor){
start = as.numeric(sfourier)*order*2
base.vec = rep(0,nrow(X)+pred_t)
base.vec[1+S*(c(1:floor((nrow(X)+pred_t)/S))-1)] = 1
for (j in 1:S){
sea[start+j,] = c(rep(0,j-1),base.vec[1:(nrow(X)+pred_t-j+1)])
}
}
W.list[[paste(c('sea',i),collapse='')]] = sea
}
}

# step 5: construct para_panel
S_n = ncol(combn(c(1:length(S.pick)),r.nomial))
para_panel = data.frame(matrix(nrow=length(p_range)*length(q_range)*S_n,ncol=7))
colnames(para_panel) = c('p','q','S1','S2','W','ic','dim')

comb = combn(c(1:length(S.pick)),r.nomial)
S_vec = c()
W_vec = c()
for (i in 1:S_n){
S_vec[i] = paste(unlist(S.pick[as.vector(comb[,i])]),collapse=' ')
W_vec[i] = paste(c(paste('sea',as.vector(comb[,i]),sep=''),'exo'),collapse=' ')
}
para_panel[,5] = rep(W_vec,length(p_range)*length(q_range))
if (sar){
para_panel[,3] = rep(S_vec,length(p_range)*length(q_range))
}else{
para_panel[,3] = NA
}
if (sma){
para_panel[,4] = rep(S_vec,length(p_range)*length(q_range))
}else{
para_panel[,4] = NA
}
para_panel[,1] = rep(p_range,rep(S_n*length(q_range),length(p_range)))
para_panel[,2] = rep(rep(q_range,rep(S_n,length(q_range))),length(p_range))

# step 6: fit and select optimal
para_panel = parallel.fit(X,W.list,para_panel,U,'AIC',multicore)

if (crit=='BC'){
dim_up = para_panel[which.min(para_panel[,6]),'dim']
for (i in 1:nrow(para_panel)){
dim = para_panel[i,'dim']
para_panel[i,'ic'] = para_panel[i,'ic'] - 2*dim + (nrow(X)^(2/3))*(sum(1/c(1:dim)))
}
temp_para_panel = para_panel[para_panel\$dim<=dim_up,]
best.row = temp_para_panel[which.min(temp_para_panel[,6]),1:5]
}else{
if (crit=='BIC'){
para_panel\$ic = para_panel\$ic - 2*para_panel\$dim + log(nrow(X))*para_panel\$dim
}
best.row = para_panel[which.min(para_panel[,6]),1:5]
}

p = best.row[1,1]
q = best.row[1,2]
if (is.na(best.row[1,3])){S1=NULL
}else{S1 = as.integer(strsplit(best.row[1,3],' ')[[1]])}
if (is.na(best.row[1,4])){S2=NULL
}else{S2 = as.integer(strsplit(best.row[1,4],' ')[[1]])}
W.temp = rbindlist(W.list,best.row[1,5])
best.m = garma.fit(X,U,p,q,S1,S2,W.temp,crit)
prediction = garma.pred(best.m,X,W.temp,pred_t)

output = list(obj=best.m,ic_panel=para_panel,prediction=prediction)
}
```
JieGroup/tsms documentation built on Sept. 15, 2020, 10:39 a.m.