# Calculation of the stationarity condition for a given model.
# This program calculate the stationary criterion for bilinear model
# with the coefficient of X(t-i)*e(t-j) equal 0 if i<j.
#
# Converted from a FORTRAN program written by
# B.K. Stensholt (1989), Statistical analysis of multivariate 
#      bilinear time series models, 
#      PhD thesis Victoria University of Manchester, UK
#
######################################################################
library(fastmatrix)   # for kronecker.prod

sum1  = array(0,c(4,4,1))
sum2  = array(0,c(4,4,1))
sum3  = array(0,c(4,4,1))
gamma = array(0,c(4,4,1))

a     = array(0,c(2,2,1))
b     = array(0,c(2,2,2,1,1)) 
sigma = array(0,c(2,2))
biga  = array(0,c(2,2))
aka   = array(0,c(4,4))
bigb  = array(0,c(2,2,2,1))
b1    = array(0,c(2,2))
b2    = array(0,c(2,2))
bkb   = array(0,c(4,4))
sbkb  = array(0,c(4,4,2,2,1))
aa    = array(0,c(2,2))
am    = array(0,c(2,2))
bl    = array(0,c(2,2))
bs    = array(0,c(2,2))
aabs  = array(0,c(2,2))
bkab  = array(0,c(4,4))
abkb  = array(0,c(4,4))
sbkab = array(0,c(4,4,2,2,1,1))
sabkb = array(0,c(4,4,2,2,1,1))
biggm = array(0,c(4,4))
rr    = array(0,c(4))
ri    = array(0,c(4))
	
# file name1 is the file of the bilinear coefficients
# write(*,*)'enter file name'

ndim = 2
ipb  = 1
iqb  = 1
ips  = 1
ns   = ndim*ips
ns2  = ns*ns
ns2q = ns2*iqb

for (l in 1:iqb){  
   for  (i in 1:ns2){
     for (j in 1:ns2){
     	sum1[i,j,l]  = 0.
	sum2[i,j,l]  = 0.
	sum3[i,j,l]  = 0.
	gamma[i,j,l] = 0.
     }
   }
}
   
biga[1,1] = 0.2
biga[1,2] = 0.3
biga[2,1] = 0
biga[2,2] = -0.5
       
aka <- kronecker.prod(biga, biga)

bigb[1,1,1,1] = 0.2
bigb[1,2,1,1] = -0.1
bigb[2,1,1,1] = 0.4
bigb[2,2,1,1] = -0.3

bigb[1,1,2,1] = 0.0
bigb[1,2,2,1] = 0.3
bigb[2,1,2,1] = -0.25
bigb[2,2,2,1] = 0.4

sigma[1,1] = 4
sigma[1,2] = 0.5
sigma[2,1] = 0.5
sigma[2,2] = 4
	
for (u1 in 1:ndim){   
  for (u2 in 1:ndim){ 
    for (i in 1:ns){  
      for (j in 1:ns){ # 52
	b1[i,j] = bigb[i,j,u1,1]
   	b2[i,j] = bigb[i,j,u2,1]
      }    
    }      
      
  bkb =  kronecker.prod(b1,b2)
  for (i in 1:ns2){   
    for (j in 1:ns2){ 
	sbkb[i,j,u1,u2,1] = bkb[i,j]*sigma[u1,u2]
 	sum1[i,j,1]       = sum1[i,j,1]+sbkb[i,j,u1,u2,1]
     } 
  }    
  
 }      # end big loop
}       

# for (l in 2:iqb){ # do big loop 60; NOTE iqb = 1
  l = 1
  for (u1 in 1:ndim){     
    for (u2 in  1:ndim){  
    
      for (i in 1:ns){    
        for (j in 1:ns){  
          bl[i,j] = bigb[i,j,u1,l]
   	  b2[i,j] = bigb[i,j,u2,l]
        }  
      }    
      
      bkb = kronecker.prod(b1,b2)
	
      for (i in 1:ns2){   
        for (j in 1:ns2){ 
	  sbkb[i,j,u1,u2,l] = bkb[i,j]*sigma[u1,u2]
  	  sum1[i,j,l]       = sum1[i,j,l]+sbkb[i,j,u1,u2,l]
        } 
      }   
      
     }    
   }      

  for (i in 1:ns){    
    for (j in 1:ns){  
   	aa[i,j] = biga[i,j]
    }   
  }     
 
  for (u1 in 1:ndim){          
     for (u2 in 1:ndim){       
#	for (s in 1:l-1){      # CHANGE for l>1
	for (s in 1:1){        #   ,,
	
         biga = am*aa
 	    for (i in 1:ns){    
	      for (j in 1:ns){  
        	aa[i,j] = am[i,j]
             }     
          }        
              
	for (i in 1:ns){    
	  for (j in 1:ns){  
	    bl[i,j] = bigb[i,j,u1,l]
  	    bs[i,j] = bigb[i,j,u2,s]
  	  }
  	}   
  	
        bs   = aabs*aa
        bkab = kronecker.prod(bl,aabs)
        abkb = kronecker.prod(aabs,bl)

	for (i in 1:ns2){  
          for (j in 1:ns2){
	    sbkab[i,j,u1,u2,l,s] = bkab[i,j]*sigma[u1,u2]
            sabkb[i,j,u1,u2,l,s] = abkb[i,j]*sigma[u1,u2]
            sum2[i,j,l] = sum2[i,j,l]+sbkab[i,j,u1,u2,l,s]
            sum3[i,j,l] = sum3[i,j,l]+sabkb[i,j,u1,u2,l,s]
          }  
       }     
       
     }    
    }     
  }       

       for (i in 1:ns2){   
         for (j in 1:ns2){ 
          gamma[i,j,l] = sum1[i,j,l]+sum2[i,j,l]+sum3[i,j,l]
         }
       }    
 
# }  # end 60 for (l in 2:iqb) NOTE iqb = 1     

for (i in 1:ns2){    
  for (j in 1:ns2){  
    gamma[i,j,1] = sum1[i,j,1]+aka[i,j]
  } 
}   

for (l in 1:iqb){      
  for (i in 1:ns2){    
    for (j in i:ns2){  
	jg = (l-l)*ns2+j
 	biggm[i,jg] = gamma[i,j,l]
    }
  }
}

# for (i in ns2+1:ns2q){  # 93, NOTE: ns2 = 4 and ns2q = 4
#   for (j in 1:ns2q){    # 93
#	biggm[i,j]=0.
#   	biggm[i,i-ns2]=1.
#   }   # 93
# }     # 93

eig   = eigen(biggm)
sqeig = sqrt(Re(eig$values)^2 + Im(eig$values)^2)  

return(sqeig)

########## END ST.R    ######################