# File: Scmi-cpi.R 12/03/2024
# This program computes the two prediction intervals:
# SCMI and CPI (blocked = MCDR; ## !! ) for a given coverage probability "alpha"
#
# De Gooijer, J.G. and Zerom, D. (2000).
#      Kernel based multi-step-ahead prediction of the U.S. short-term interest
#      rate. Journal of Forecasting, 19(4), 335-353.      
##############################################################################
# constructing the product kernel using Gaussian distribution 

ker = function(x){(exp(-0.5*x))*((2*pi)^(-k/2))}


### This supporting function is the cumulative density ########
cum = function(y1,h1=hop){   # 1)

sum1 = 0     
sum2 = 0
t1   = 1

while ( t1 <= n0){     
    vect1 = (z[n, ]-x0[t1,])/h1       
    k2    = t(vect1) %*% (vect1) 
    sum1  = sum1 + ker(k2) 
    if (any(yy[t1,] <= y1)){
      vect2 = (z[n,]-x0[t1,])/h1    
      k1    = t(vect2) %*% (vect2)
      sum2  = sum2 + ker(k1)
    }
    t1    = t1+1     
}

if (sum1 == 0){
  sum1 = 1.0e-300
}

return(sum2/sum1)
}    # end 1) 


###########  MAIN  #################################
tbill9 <- read.table("C:/ali-third/Computer Code/mean-mode-median-Gauss/tbill9.txt", quote="\"", comment.char="")
data   = tbill9
datats = ts(data,freq=1)+10

Miny = min(tbill9)
Maxy = max(tbill9)
yy   = data # tbill9

n1     = nrow(data)
data1  = datats  
      
k = 6      # The Markov coefficient  CHANGE

## Construct lagged values
x1 = data1[6:n1,]     
x2 = data1[5:(n1-1),]     
x3 = data1[4:(n1-2),]     
x4 = data1[3:(n1-3),]     
x5 = data1[2:(n1-4),]     
x6 = data1[1:(n1-5),] 

## Alternative for the above x1, ... ,x6
## data1 = {}
##  j = 0
##  while ( j<= (k-1) ){
##    data2 = data[k-j:n1-j,]
##    data1 = cbind(data1,data2)
##    j = j+1
##  }


data = cbind(x1,x2,x3,x4,x5,x6)     
z    = data     
        
zn = {}     
n1 = {}  
## mcdr2 ={}

width11 = matrix(0,1,1)    
width22 = matrix(0,1,1)     

zo1     = {}      
cpo3    = {}     
count0  = {}     
count10 = {} 

M_scm2 = matrix(0,5,2)  # shortest conditional modal interval
M_cpo2 = matrix(0,5,2)  # shortest percentile CI

mx  = 5    # maximum number of forecasts
prd = 100  # 254  # number of replications, CHANGE 
m   = 1     

while (m <= mx){    # 1) start BIG loop  

## ii1 = {}

  r1   = {}     
  cpo  = {}     
  scm1 = {} 
## mcdr1 = {}

  tru    = {}     
  zo     = {}     
  count1 = 0     
  count  = 0  
  
  cpo2  = {}     
  scm2  = {}  
## mcdr2 = {}

  n   = 1670   # (just an example ) the starting time index 
  nob = n + prd - 1
    
  while (n <= nob ){ #  1923){  # 2) start big loop
  ii = n - 1670 + 1
  
  mn = cbind(ii,m,n)
  print("ii,  m and n = ")
  print (mn)

     alpha = 0.9  # Coverage probability. Just an example. CHANGE 
      
     NK  = 20     
     n0  = n-m     
     y   = z[m+1:n,1]     
     x0  = z[1:n0,]     
     std = sd(data1[1:n])     
     hop = std*(n)^(-1/(k+4)) 
#
#    Here we calculate the cumulative density
#    in a grid of 256 values 
     F  = rep(0,256)     
     y2 = rep(0,256)     
     a  = Miny-3*hop     
     b  = Maxy+3*hop     
     
     step = (b-a)/256     
     i = 1     
     while (i <= 256){ 
        y2[i] = a+(i-1)*step     
        F[i]  = cum(y2[i],hop)     
        i = i+1     
      }  
      
#     We use from now on the matrix "matt"       
      matt = cbind(y2,F)     
      r    = nrow(matt)     

# (1) Computation of the percentile CI (cpo) 
     j = 1     
     while (j <= r ){     
       if ( matt[j,2] >= (0.5-(alpha/2)) ) {    
         low = matt[j,1]     
         break     
        }
        j = j+1     
     }

     # The first part of the loop works as follows. As
     # soon as that value is obtained, it is stopped 

     up1 = {}     
     o   = j+1     
     while ( o <=r ){     
       if ( matt[o,2] >= (0.5 + (alpha/2)) ){     
          upper1 = cbind(matt[o,1],matt[o,2])     
          up1    = rbind(up1,upper1)     
       }
       o = o+1     
      }

       # This looks for ties so that it
       # only returns the last value.  Also for
       # very small
       # differences in probability values, "i" extend
       # the value further to the right 

       u       = nrow(up1)     
       upperbo = matrix(0,1,1)     

        if (u != 1){                                                                                                                                                     
          tol = 1.0e-6     
          f   = 2     
          while (f <= u){     
             if ( abs(up1[f,2] - up1[f-1,2]) > tol){     
               upperbo = up1[f-1,1]     
               break     
             }     
             f = f+1     
           }
           if (upperbo == 0){     
              upperbo = up1[u,1]     
            }

         }

       if ( u== 1){     
           upperbo = up1[1,1]     
       }

       inter = cbind(low,upperbo) # percentile confidence interval   
       cpo1  = rbind(cpo1,inter)
       
########################################
# (2) Computation of the shortest conditional modal (scm) interval 

   minval   = 2000     
   cover    = 0     
   ydif     = 0     
   gap      = 0     
   scm      = matrix(0,1,2)     
   gg       = 0     
   FF       = matt     
   scmlower = 0     
   scmupper = 0     
   toler    = 1.0e-5     

   i = 1     
   while (i <= r ){            # 4) start  
     if  ( FF[i,2] <= 0.3 ){   # 5) start   
       j = i+1     
       while (j <= r ){        # 6)  start
         gap = FF[j,2]-FF[i,2]     
         if ( gap >= alpha ){  # 7) start  
           ydif = FF[j,1]-FF[i,1]     
           if (( ydif <= minval ) | ((ydif==minval) & (gap>cover))){      
             minval   = ydif     
             cover    = gap     
             scmlower = FF[i,1]     
             scmupper = FF[j,1]     
           }   # end if
           break     
         }       # end 7)
         j = j+1     
       }         # end 6)    

     } else if ( FF[i,2] < 0.075 ){    
       i = i+1     
     }           # end 5)

       else {     
      break     
   }  # end   

   i = i+1     
 }               # end 4)

 scm  = cbind(scmlower,scmupper) #shortest cond. modal (scm) interval    
 scm1 = rbind(scm1|scm)

################# BLOCKED BY THE NOTATION "##" ###################### 
# (3) Computation of the maximum conditional density region. 
# This is computationally heavy. It requires at least 
# 5 minutes on its own. 

## mcdr = matrix(0,1,4)
## tmp2 = 0
## scl1 = 0     
## scu1 = 0     
## scl2 = 0     
## scu2 = 0     

## minval = 1000     
## P  = 0     
## i1 = 0     
## i2 = 0     

## i = 1     
## while( i <= (r-1) ){  # start SUPER BIG loop    
##    if ( FF[i,2] <= alpha1){     
##         i = i+1     
##    }     
##    end if {     
##
##    if (FF[i,2] > (1-alpha)){     
##       break     
##    }      
##
##   j = i+1     
##   while (j <= (r-1) ){    # start BIG loop     
##     e = j     
##     while ( e <= (r-1) )  # start big loop      
##
##      l = e+1     
##      while (l <= (r-1) ){       # 0)   
##        gap = FF[j,2]-FF[i,2]+FF[l,2]-FF[e,2]     
##        if (gap >= alpha){       # 1)  
##          tmp2 = FF[j,1]-FF[i,1]+FF[l,1]-FF[e,1]     
##          if ( tmp2 < minval){   # 2)   
##            scl1   = FF[i,1]     
##            scu1   = FF[j,1]     
##            scl2   = FF[e,1]     
##            scu2   = FF[l,1]     
##            minval = tmp2     
##            p  = gap     
##            i1 = FF[j,2]-FF[i,2]     
##            i2 = FF[l,2]-FF[e,2]     
##
##          } else if ((tmp2==minval) & (gap>p)){     
##            scl1   = FF[i,1]     
##            scu1   = FF[j,1]     
##            scl2   = FF[e,1]     
##            scu2   = FF[l,1]     
##            minval = tmp2     
##            p  = gap     
##            i1 = FF[j,2]-FF[i,2]     
##            i2 = FF[l,2]-FF[e,2]     
##        } # end if 2)
##       break     
##       }  # end if 1)
##       
##       l = l+1     
##       }  # end while 0)     
##
##     e = e+1     
##     }  # end big loop     
##   j = j+1     
##   }    # end BIG loop     
##
##   # else if (FF[i,2] < 0.025){     
##   #   i = i+1     
##   #   } else {    
##   #   break     
##   # }   
##   
##  i = i+1     
##  } # end SUPER BIG loop      
##
##
##  mcdr = cbind(scl1,scu1,scl2,scu2)     
##  ii   = cbind(i1,i2)     
##
##   mcdr1 = rbind(mcdr1,mcdr)     
##   ii1   = rbind(ii1,ii)      
 
 
###############  end MCDR #############################

   if ( (yy[n+m,] <= scm[1,2]) & (yy[n+m,] >= scm[1,1]) ){     
     coun  = 1     
     count = count+coun     
     scm1  = rbind(scm1,scm)     
   }    

   if ( (yy[n+m,] <= inter[1,2]) & (yy[n+m,] >= inter[1,1]) ){     
     coun1  = 1     
     count1 = count1+coun1     
     cpo    = rbind(cpo,inter)     
    }

    scm2 = rbind(scm2,scm)  # shortest cond. modal (scm) interval      
    cpo2 = rbind(cpo2,inter)     

    n = n+1     
  }  # end 2) big loop  
   
  print("mean scm2 = ")
  M_scm2[m,] = colMeans(scm2)  
  print(colMeans(scm2))
  
  print("mean cpo2 = ")
  M_cpo2[m,] = colMeans(cpo2)  
  print(colMeans(cpo2))
   
###################################################

   # The average width for only those cases where
   # the true value is incorporated 

    widthw1 = sum(abs(cpo[,2]  - cpo[,1]))*(1/(count1))     
    widthw2 = sum(abs(scm1[,2] - scm1[,1]))*(1/(count))     

   # Counts the number of instances where the true value
   # is incorporated 

   count0  = rbind(count0,count)     
   count10 = rbind(count10,count1)     

   width11 = widthw1 
   width11 = rbind(width11,widthw1)  # cpo, percentile CI
   width22 = rbind(width22,widthw2)  # scm, shortest cond modal   

#   scm3[m,] = scm2     
#   cpo3 = cbind(cpo3,cpo2)     

  m = m+1     
} # 1) end BIG loop  

#########################################################
#########################################################

print("n = ")
print (n)
print("prd = ")
print(prd)

# M_scm2, for 100 replications, AR(6) model
#           [,1]      [,2]
#[1,] -0.1132721 0.1223581
#[2,] -0.1132721 0.1223581
#[3,] -0.1134453 0.1221848
#[4,] -0.1137930 0.1218372
#[5,] -0.1136199 0.1218370
# M_cpo2
#           [,1]      [,2]
#[1,] -0.1061568 0.1249663
#[2,] -0.1061568 0.1249663
#[3,] -0.1065032 0.1249663
#[4,] -0.1066767 0.1249663
#[5,] -0.1065036 0.1249663