#  This code gastma.R   29/02/2024
#  Tests linearity of an MA(1) against GASTMA(1) model.
#  MA(1) parameter of the model estimated using the Hannan-Rissanen
#  technique. This version assumes constant unknown and estimates it.
#  Reference:Brannas, K., De Gooijer, J.G., and Terasvirta, T. (1998),
#             "Testing linearity against nonlinear moving average models",
#    Communications in Statistics -  Theory and Methods, 27(8), 2025--2035
#  ************************************************************************ 

# 1)   Start the huge loop 
lambda = 1     
khuge  = 0     

while (khuge<=1){  # 1) start huge loop   

#  Give the number of observations 
if (khuge==1){
  norig = 51
  ny    = 5     
} else{
  norig = 201    
  ny    = 5     
}

# 2)  Start the big loop
if (khuge==1){     
   kbig = 0     
}  else {     
   kbig = 0     
}

theta = 0.8     

while (kbig<=2){    # 2) Start the big loop

  if (kbig==1){      
     theta = -0.8      
   } else if (kbig==2){      
     theta = 0      
   } else {      
     theta = 0.8      
 }  # end if  

#  Generate the random numbers 
# rndseed = 430308     
knt     = 1     
kntr    = 1     
pvknt   = matrix(0,9,1)     
kntmax  = 1000     
rmatrix = matrix(0,norig,kntmax)     

# 3)  Start the loop 

while (knt<=kntmax){     
# Generate the random normal numbers 
u              = rnorm(norig,0,1)     
rmatrix[ ,knt] = u     

ulag1 = u[1:(norig-1)]
u     = u[2:norig]     
yxx   = matrix(0,(norig-1),1)     
yxxx  = cbind(yxx,ulag1)     

uplus1 = apply(yxxx,2,max)  # find maximum

#  Create y 
bknt = 1     
if (kbig==1){     
  alpha = 0     
  } else if (kbig==2){     
  alpha = -0.8     
  } else {    
  alpha = -1.6     
}

while (bknt <= 9){  # start 4)
   #  y=u+theta*ulag1+alpha*uplus1; 
   # Y_t = e_t + theta*e_{t-1} + alpha * e_{t-1}*( 1/exp(-e_{t-1}) )
   y = u+theta*ulag1+alpha*(1./(1+exp(-lambda*ulag1)))*ulag1 
   y<-y[!is.na(y)]  # remove na in vector
   n = nrow(as.matrix(y))  
   # x = matrix(0,(n-ny-1),1)   
   AR[bknt] = alpha
   TH[bknt] = theta

   # Estimate the MA(1) model using Hannan-Rissanen with an AR(ny=5) model
   #   (1) Create lags of y
   j    = 1  
   ylag = matrix(0,(n-ny),ny)    

   while (j <= ny){
     bbb = y[(ny-j+1):(n-j)]  # j=1, bbb[5-1+1:n]=y[5:500] 
      #  j=2, bbb=y[5-2+1:200]=y[4:500]
      ylag[,j] = bbb   
      j = j+1     
   } 
   n1 = nrow(as.matrix(bbb))
   
   yy = y[ny+1:n]    
   y  <- yy[!is.na(yy)]  # remove NAs

   ry   = nrow(as.matrix(y))     
   yks  = matrix(1,ry,1)     
   ylag = cbind(yks, ylag)

   # Estimate an AR(ny) model 
   # library(matlib)
   minv = solve(t(ylag) %*% ylag)
   uha  = y - ylag %*% minv %*%t(ylag) %*% y     

   uha1 = uha[1:n1-1]     
   y    = y[2:n1]   
   ry   = nrow(as.matrix(y))     
   yks  = matrix(1,ry,1)     

   uha1 = cbind(yks,uha1)     

   # Estimate a MA(1) model 
   bhat = solve(t(uha1) %*% uha1) %*% t(uha1) %*% y     

   if (abs(bhat[2,1])>1){     
     bhat[2,1]=1/bhat[2,1]     
   }
  
   uhat = y-uha1 %*% bhat     
   ruh  = nrow(as.matrix(uhat))

   #  Store bhat[2,1]    
   bha = bhat[2,1]     
   if (bknt==9){     
     if (knt==1){     
       b2cum = bha     
     } else {     
     b2cum = rbind(b2cum,bha)     
    }     
   }   # end if 

   # Create a series of squared errors 
   ypos = uhat^2     
 
   #  Trim the series 
   bhat = bhat[2,1]  
 
   if (abs(bhat)<0.1){     
    jlag = 1     
    u1lags = uhat[1:(ry-1)]     
    y1lags = ypos[1:(ry-1)]     
   } 
  
   if(abs(bhat^2)<0.1){
    jlag = 2     
    u1lags = uhat[2:(ry-1)]+bhat*uhat[1:(ry-2)]    
    y1lags = ypos[2:(ry-1)]+bhat*ypos[1:(ry-2)]     
   }  
  
   if (abs(bhat^3)<0.1){     
    jlag = 3     
    u1lags = uhat[3:(ry-1)]+bhat*uhat[2:(ry-2)]+bhat*bhat*uhat[1:(ry-3)]     
    y1lags = ypos[3:(ry-1)]+bhat*ypos[2:(ry-2)]+bhat*bhat*ypos[1:(ry-3)]     
   } 
  
   if (abs(bhat^4)<0.1){     
    jlag = 4     
    u1lags = uhat[4:(ry-1)]+bhat*uhat[3:(ry-2)]+bhat*bhat*uhat[2:(ry-3)]+(bhat^3)*uhat[1:(ry-4)]     
    y1lags = ypos[4:(ry-1)]+bhat*ypos[3:(ry-2)]+bhat*bhat*ypos[2:(ry-3)]+(bhat^3)*ypos[1:(ry-4)]     
   } 
   
   if (abs(bhat^5)<0.1){       
    jlag = 5     
    u1lags = uhat[5:(ry-1)]+bhat*uhat[4:(ry-2)]+bhat*bhat*uhat[3:(ry-3)]+(bhat^3)*uhat[2:(ry-4)]+(bhat^4)*uhat[1:(ry-5)]     
    y1lags = ypos[5:(ry-1)]+bhat*ypos[4:(ry-2)]+bhat*bhat*ypos[3:(ry-3)]+(bhat^3)*ypos[2:(ry-4)]+(bhat^4)*ypos[1:(ry-5)]     
   }     

   y    = y[(jlag+1):ry]     
   ryy  = nrow(as.matrix(y)) 
   uhat = uhat[(jlag+1):ruh]     
   yks  = matrix(1,ryy,1)     

   #  Run the regression 
   u11 = cbind(yks,u1lags)     
   z1  = cbind(yks, u1lags, y1lags)  

   # er=y-u11*inv(u11'u11)*u11'y     
   er   = uhat - u11 %*% solve(t(u11) %*% u11) %*% t(u11) %*% uhat    
   eur  = er - z1 %*% solve(t(z1) %*% z1) %*%t(z1) %*% er     
   sgma = t(eur)%*%eur/(ryy-3)     

  #  Run the F test 
  F  = (t(er) %*% er - t(eur)%*%eur)/sgma     
  pv = 1 - pf(F,1,(ryy-3))   
  
  # pv= 1 - pf(F,1,rows(y)-3+jlag)     

  if (pv<0.05){       # 0.5 = nominal significance level
    pvknt[bknt,1] = pvknt[bknt,1]+1     
  }

  #  Next alfa 
  bknt  = bknt+1
  alpha = alpha+0.2  
  
}  # end 4) loop (bknt <= 9) 
################################

  knt  = knt+1    
  kntr = kntr+1     
  if (kntr>100){
    kntr = 1
  }

}  # 3) end loop knt<kntmax  

cat("File:gastma6.R", "\n")  
cat("Tests linearity of an MA(1) against GASTMA(1)", "\n") 

cat("\n")
if (kbig==1){     
    print("The number of observations is: norig-1")
    print("Theta = -0.8, Alpha = 0, 0.2,...,1.6 (9 values)")
  } else if (kbig==2){     
    print("The number of observations is: norig-1 ")
    print("Theta = 0, Alpha = -0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0,6,0.8 (9 values) ")     
  } else {     
    print("The number of observations is: norig-1 ")
    print("Theta = 0.8, Alpha = -1.6,-1.4,-1.2,-1.0,-0.8,-0.6,-0.4,-0.2,0 (9 values) ")     
}

cat("\n")
print("Number of trials is: knt-1, and the p-values are:")
print(knt-1)
pval = pvknt/(knt-1)     
print(pval) 

p_tot = cbind(AR,pval)

# Print the mean and std of bhat[2,1] and save b2cum 
b2mean = mean(b2cum)     
b2std  = sd(b2cum)  

cat("\n")
print("Mean and standard deviation of b[2,1] are: ")     
print(b2mean)
print(b2std)

kbig = kbig+1     

}  # end the big loop 2)    

# End of the huge loop 
khuge = khuge+1     

}  # end 1) huge loop