#source this AFTER sourcing wtdchisqfun.r; needed on penultimate line below

multisteptest=function(sarimamod,lead){
#Returns the significance level of a statistic for testing whether 
#a significant reduction in multistep forecast error variance, at the
#specified lead time, may be achieved by adjusting the parameters of
#the specified arima model which has been fitted by maximum likelihood
#using the arima function.

#An internal function for convolution is required
conv=function(v)
{
n=length(v)-1
u=rep(0,(n+1))
for (i in 0:n)
{
s=c(0)
for (j in 0:(n-i))
{
s=s+v[j+1]*v[j+1+i]
}
u[i+1]=s
}
u
}
#The required structures are extracted from the arima model 
pdq=sarimamod$arma
p=pdq[1]
q=pdq[2]
P=pdq[3]
Q=pdq[4]
s=pdq[5]
d=pdq[6]
D=pdq[7]

a=sarimamod$coef
if (p>0) {phi=a[1:p];phim=c(1,-phi)} else {phim=1}
if (q>0) {theta=-a[(1+p):(q+p)];thetam=c(1,-theta)} else {thetam=1}
if (P>0) {Phi=a[(1+p+q):(P+p+q)];Phim=c(1,-Phi)} else {Phim=1}
if (Q>0) {Theta=-a[(1+p+q+P):(Q+p+q+P)];Thetam=c(1,-Theta)} else {Thetam=1}
#phi;theta;Phi;Theta
cphi=conv(phim)
cphi
ctheta=conv(thetam)
ctheta
cPhi=conv(Phim)
cTheta=conv(Thetam)
maxpq=max(p,q)
maxPQ=max(P,Q)
stt=1+maxpq+maxPQ*s+d+D*s
nreg=p+q+P+Q
xres=tsmod$residuals  #the test uses the model residuals
#omit early orthogonal residuals that are not stationary
#or missing values due to start up of kalman filter 
xres=xres[stt:length(xres)] 
sigsq=tsmod$sigma2
nres=length(xres)
halfn=floor(nres/2)+1
hnp1=halfn+1
find=(1:nres)
ones=rep(1,nres)
freq=(find-1)/nres
#Construct frequency domain response and regressors
Y=(abs(fft(xres))^2)/(nres*sigsq) -1
if (maxpq>0){
cvec=array(0,c(nres,maxpq))
for (k in (1:maxpq)) {cvec[,k]=2*cos(2*pi*freq*k)}  
}
if (maxPQ>0){
cvecs=array(0,c(nres,maxPQ))
for (k in (1:maxPQ)) {cvecs[,k]=2*cos(2*pi*freq*k*s)}  
}
if (p>0) {
  Xp=cbind(cvec[,(1:p)])
  M=cbind(ones,Xp)
  phiden=M%*%cphi
  for (k in (1:p)) {Xp[,k]=Xp[,k]/phiden}
}
if (q>0) {
  Xq=cbind(cvec[,(1:q)])
  M=cbind(ones,Xq)
  thetaden=M%*%ctheta
  for (k in (1:q)) {Xq[,k]=Xq[,k]/thetaden}
}
if (P>0) {
  XP=cbind(cvecs[,(1:P)])
  M=cbind(ones,XP)
  Phiden=M%*%cPhi
  for (k in (1:P)) {XP[,k]=XP[,k]/Phiden}
}
if (Q>0) {
  XQ=cbind(cvecs[,(1:Q)])
  M=cbind(ones,XQ)
  Thetaden=M%*%cTheta
  for (k in (1:Q)) {XQ[,k]=XQ[,k]/Thetaden}
}
Xreg=array(0,c(nres,nreg))
Z=Xreg
if (p>0) {Xreg[,(1:p)]=Xp}
if (q>0) {Xreg[,(p+(1:q))]=Xq}
if (P>0) {Xreg[,(p+q+(1:P))]=XP}
if (Q>0) {Xreg[,(p+q+P+(1:Q))]=XQ}
#Orthogonal regressors to improve conditioning without affecting the test
Xsvd=svd(Xreg)
X=Xsvd$u
#Construct model psi-weights
impul=rep(0,nres)
impul[1]=1
psiwt=impul
psiwts=impul
for (t in (1:lead))
 {sum=impul[t]
  if (q>0){for (k in (1:q)){if ((t-k)>0) {sum=sum-theta[k]*impul[t-k]}}}
  if (p>0){for (k in (1:p)){if ((t-k)>0) {sum=sum+phi[k]*psiwt[t-k]}}}
  psiwt[t]=sum}
for (t in (1:(lead+1)))
 {sum=psiwt[t]
  if (Q>0){for (k in (1:Q)){if ((t-k*s)>0) {sum=sum-Theta[k]*psiwt[t-k*s]}}}
  if (P>0){for (k in (1:P)){if ((t-k*s)>0) {sum=sum+Phi[k]*psiwts[t-k*s]}}}
  psiwts[t]=sum}    
if (d>0){
  for (i in (1:d) ){
    for (t in (1:(lead+1))) {if ((t-1)>0) {psiwts[t]=psiwts[t]+psiwts[t-1]}}}}
if (D>0){
  for (i in (1:D) ){
    for (t in (1:(lead+1))) {if ((t-s)>0) {psiwts[t]=psiwts[t]+psiwts[t-s]}}}}
#plot(psiwts)  
psitran=fft(psiwts)
lagmask=rep(1,nres)
lagmask[(1:lead)]=rep(0,lead)
lagmask[(hnp1:nres)]=rep(0,(nres-hnp1+1))
#Construct modified regressors
for (j in (1:nreg)){
 zvec=psitran*X[,j]
 zvectran=fft(zvec,inverse=TRUE)*lagmask
 zvec=fft(zvectran)/nres
 Z[,j]=2*Re(zvec*Conj(psitran))
}
#Construct test statistic and distribution weights by well conditioned method
Zsvd=svd(Z)
Zu=Zsvd$u
Zd=Zsvd$d
Zv=Zsvd$v
Hi=diag(Zd)%*%(solve((t(Zu)%*%X)%*%Zv))
g=t(Zu)%*%Y
test=as.numeric((t(g)%*%Hi%*%g)*sigsq/2)
q=test*2/nres
A=Zu-X%*%(t(X)%*%Zu)
Asvd=svd(A)
RT=Asvd$v%*%diag(Asvd$d)
S=t(RT)%*%Hi%*%RT
D=(eigen(S)$values)*sigsq
probH=Pwtdchisq(test,D)
list(statistic=test,weights=D,pvalue=probH)
# end of function multisteptest
}