### This code applies the Zhang and Xiang (2008) methodology to obtain tail probabilities (Equations (4)-(6) in the paper)
###
###
###
###
###
rm(list = ls())
### change some basic options in R
options(java.parameters = "-Xmx8000m")
options( digits = 7, scipen = 8 )
options( prompt ="R> ", width = 64)
options( show.signif.stars = FALSE )
### load necessary libraries
library( "animation" )
library( "dynlm" )
library( "fOptions" )            # to calculate the nominal price of the RR
library( "hydroTSM" )
library( "lmtest" )
library( "NMOF" )
library("openxlsx")             # no more problems with Java
library( "pracma" )
library( "sandwich" )
library( "tseries" )            # kpss.test
library( "xtable" )
library( "xts" )
### parameters
linestoskip <- 0
### change working directory
data.file.dir <- "//XXX/Eigene Dateien/Forschung/Daniel Nathan"
setwd( data.file.dir )
### load the required data files
### read data file (i.e. the butterfly spread, risk reversal and the atm TS)
  data.file.dir <- "//XXX/Eigene Dateien/data/implicit TZ/TZ ILS.xlsx"
  data <- read.xlsx( data.file.dir, "data", startRow = 2, detectDates = TRUE )
  for( i in 2:ncol(data) ){
    data[,i] <- xts( data[,i], order.by = data[,1] )
  }
  data <- data[,-1]
  ### remove 2m-data
    data <- subset( data, select = -c( 2, 8, 14, 20, 26, 37, 43, 50 ) )
  ### choose only necessary data
    data <- subset( data, select = c(1:26, 31:40, 42:46 ) )
  ### 
    data <- na.omit( data )
```



```{r}
IV_smile <- function( gamma0, d, par ){
  gamma1 <- par[1]
  gamma2 <- par[2]
  IV <- gamma0 * ( 1 + gamma1 * d + gamma2 * d^2 )
  return(IV)
}
meansq_price_error <- function( gamma0, data, par ){
  gamma0 <- gamma0
  gamma1 <- par[1]
  gamma2 <- par[2]
  # IV_smile <- gamma0 * ( 1 + gamma1 * d + gamma2 * d^2 )
  z <- 0.2 * ( ( data[2] - IV_smile( gamma0, d = data[2],  par ) )^2 + ( data[3] - IV_smile( gamma0, data[4], par ) )^2 + ( data[5] - IV_smile( gamma0, data[6], par ) )^2 + ( data[7] - IV_smile( gamma0, data[8], par ) )^2 + ( data[9] - IV_smile( gamma0, data[10], par ) )^2 ) 
  return( z )
}

IV <- function( gamma, d ){
  gamma0 <- gamma[,1]
  gamma1 <- gamma[,2]
  gamma2 <- gamma[,3]
  IV <- gamma0 * ( 1 + gamma1 * d + gamma2 * d^2 )
  return(IV)
}
### tail probs
A <- matrix( NA, nrow = dim(data)[1], ncol = 30 )
colnames( A ) <- c( "Drop1M2", "Drop3M5", "Drop6M8", "Drop9M10", "Drop12M10", "Drop1M3", "Drop3M6", "Drop6M9", "Drop9M11", "Drop12M12", "Drop1M4", "Drop3M7", "Drop6M10", "Drop9M12", "Drop12M14", 
                    "Rise1M2", "Rise3M5", "Rise6M8", "Rise9M10", "Rise12M10", "Rise1M3", "Rise3M6", "Rise6M9", "Rise9M11", "Rise12M12", "Rise1M4", "Rise3M7", "Rise6M10", "Rise9M12", "Rise12M14" )


### 1 month
  tt <- 1 / 12
  F <- data[,"FORW1M"]
  K_A <- F * exp( 0.5 * ( data[,"ATM1m" ]/100 )^2 * tt )
  IV_A <- data[,"ATM1m" ]
  ### 10-Delta options
    delta <- 0.1
    alpha <- -qnorm( delta * exp( log( 1 + data[,"USLIBOR1m" ] / 100 ) * tt ) )
    IV_C <- data[,"BF101m" ] + data[,"ATM1m" ] + 0.5 * data[,"RR101m" ]
    IV_P <- data[,"BF101m" ] + data[,"ATM1m" ] - 0.5 * data[,"RR101m" ]
    K_10C <- F * exp( 0.5 * ( IV_C/100 )^2 * tt + alpha *  IV_C/100 * sqrt( tt ) )
    K_10P <- F * exp( 0.5 * ( IV_P/100 )^2 * tt - alpha *  IV_P/100 * sqrt( tt ) )
    IV_10C <- IV_C
    IV_10P <- IV_P   
  ### 25-Delta options
    delta <- 0.25
    alpha <- -qnorm( delta * exp( log( 1 + data[,"USLIBOR1m" ] / 100 ) * tt ) )
    IV_C <- data[,"BF251m" ] + data[,"ATM1m" ] + 0.5 * data[,"RR251m" ]
    IV_P <- data[,"BF251m" ] + data[,"ATM1m" ] - 0.5 * data[,"RR251m" ]
    K_25C <- F * exp( 0.5 * ( IV_C/100 )^2 * tt + alpha *  IV_C/100 * sqrt( tt ) )
    K_25P <- F * exp( 0.5 * ( IV_P/100 )^2 * tt - alpha *  IV_P/100 * sqrt( tt ) )  
    IV_25C <- IV_C
    IV_25P <- IV_P  
  ### gsi
    sig <- 0.2 * ( IV_A + IV_10C + IV_10P + IV_25C + IV_25P ) 
    gsi_A <- log( K_A/F ) / ( sig/100 * sqrt( tt ) )
    gsi_10C <- log( K_10C/F ) / ( sig/100 * sqrt( tt ) )    
    gsi_10P <- log( K_10P/F ) / ( sig/100 * sqrt( tt ) ) 
    gsi_25C <- log( K_25C/F ) / ( sig/100 * sqrt( tt ) )    
    gsi_25P <- log( K_25P/F ) / ( sig/100 * sqrt( tt ) ) 
  ### matrix with the data
    M <- matrix( NA, ncol = 11, nrow = length( IV_A ) )
    M[,1] <- IV_A/100
    M[,2] <- gsi_A
    M[,3] <- IV_10C/100
    M[,4] <- gsi_10C
    M[,5] <- IV_10P/100
    M[,6] <- gsi_10P
    M[,7] <- IV_25C/100
    M[,8] <- gsi_25C
    M[,9] <- IV_25P/100
    M[,10] <- gsi_25P
    M[,11] <- as.Date( index( IV_A ) )
    colnames( M ) <- c( "IV_A", "gsi_A", "IV_10C", "gsi_10C", "IV_10P", "gsi_10P", "IV_25C", "gsi_25C", "IV_25P", "gsi_25P", "dates" )
    M <- na.omit(M)
    Par <- matrix( NA, nrow = dim( M )[1], ncol = 3 )
    Par[,1] <- M[,1]
    for( i in 1:dim(M)[1] ){
      Par[i,2:3] <- stats::optim( par = c( 0.1, 0.05 ),
                                  fn = meansq_price_error, data = M[i,], gamma0 = M[i,1],
                                  method = "L-BFGS-B" )$par #,
      #control = list(trace = 5) )
    }
    ### drop 1
      S <- 0.975 * data[, "FX" ] 
      sig_bar <- window( sig, as.Date( M[,11] ) )/100
      memo <- window( log( S/data[, "FORW1M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW1M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,1] <- F              # NULL < 0
    ### drop 2
      S <- 0.97 * data[, "FX" ]
      memo <- window( log( S/data[, "FORW1M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW1M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,6] <- F              # NULL < 0
    ### drop 3
      S <- 0.96 * data[, "FX" ]
      memo <- window( log( S/data[, "FORW1M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW1M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,11] <- F             # 21 < 0
    ### rise 1
      S <- 1.025 * data[, "FX" ]
      memo <- window( log( S/data[, "FORW1M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW1M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,16] <- 1 - F         # NULL < 0
    ### rise 2
      S <- 1.03 * data[, "FX" ] 
      memo <- window( log( S/data[, "FORW1M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW1M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,21] <- 1 - F         # NULL < 0
    ### rise 3
      S <- 1.04 * data[, "FX" ] 
      memo <- window( log( S/data[, "FORW1M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW1M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,26] <- 1 - F         # 60 < 0 
### 3 months
      tt <- 3 / 12
      F <- data[,"FORW3M"]
      K_A <- F * exp( 0.5 * ( data[,"ATM3m" ]/100 )^2 * tt )
      IV_A <- data[,"ATM3m" ]
  ### 10-Delta options
    delta <- 0.1
    alpha <- -qnorm( delta * exp( log( 1 + data[,"USLIBOR3m" ] / 100 ) * tt ) )
    IV_C <- data[,"BF103m" ] + data[,"ATM3m" ] + 0.5 * data[,"RR103m" ]
    IV_P <- data[,"BF103m" ] + data[,"ATM3m" ] - 0.5 * data[,"RR103m" ]
    K_10C <- F * exp( 0.5 * ( IV_C/100 )^2 * tt + alpha *  IV_C/100 * sqrt( tt ) )
    K_10P <- F * exp( 0.5 * ( IV_P/100 )^2 * tt - alpha *  IV_P/100 * sqrt( tt ) )
    IV_10C <- IV_C
    IV_10P <- IV_P   
  ### 25-Delta options
      delta <- 0.25
      alpha <- -qnorm( delta * exp( log( 1 + data[,"USLIBOR3m" ] / 100 ) * tt ) )
      IV_C <- data[,"BF253m" ] + data[,"ATM3m" ] + 0.5 * data[,"RR253m" ]
      IV_P <- data[,"BF253m" ] + data[,"ATM3m" ] - 0.5 * data[,"RR253m" ]
      K_25C <- F * exp( 0.5 * ( IV_C/100 )^2 * tt + alpha *  IV_C/100 * sqrt( tt ) )
      K_25P <- F * exp( 0.5 * ( IV_P/100 )^2 * tt - alpha *  IV_P/100 * sqrt( tt ) )  
      IV_25C <- IV_C
      IV_25P <- IV_P  
  ### gsi
    sig <- 0.2 * ( IV_A + IV_10C + IV_10P + IV_25C + IV_25P ) 
    gsi_A <- log( K_A/F ) / ( sig/100 * sqrt( tt ) )
    gsi_10C <- log( K_10C/F ) / ( sig/100 * sqrt( tt ) )    
    gsi_10P <- log( K_10P/F ) / ( sig/100 * sqrt( tt ) ) 
    gsi_25C <- log( K_25C/F ) / ( sig/100 * sqrt( tt ) )    
    gsi_25P <- log( K_25P/F ) / ( sig/100 * sqrt( tt ) ) 
  ### matrix with the data
    M <- matrix( NA, ncol = 11, nrow = length( IV_A ) )
    M[,1] <- IV_A/100
    M[,2] <- gsi_A
    M[,3] <- IV_10C/100
    M[,4] <- gsi_10C
    M[,5] <- IV_10P/100
    M[,6] <- gsi_10P
    M[,7] <- IV_25C/100
    M[,8] <- gsi_25C
    M[,9] <- IV_25P/100
    M[,10] <- gsi_25P
    M[,11] <- as.Date( index( IV_A ) )
    colnames( M ) <- c( "IV_A", "gsi_A", "IV_10C", "gsi_10C", "IV_10P", "gsi_10P", "IV_25C", "gsi_25C", "IV_25P", "gsi_25P", "dates" )
    M <- na.omit(M)
    Par <- matrix( NA, nrow = dim( M )[1], ncol = 3 )
    Par[,1] <- M[,1]
    for( i in 1:dim(M)[1] ){
      Par[i,2:3] <- stats::optim( par = c( 0.1, 0.05 ),
                                  fn = meansq_price_error, data = M[i,], gamma0 = M[i,1],
                                  method = "L-BFGS-B" )$par #,
      #control = list(trace = 5) )
    }
    ### drop 1
      S <- 0.95 * data[, "FX" ] 
      sig_bar <- window( sig, as.Date( M[,11] ) )/100
      memo <- window( log( S/data[, "FORW3M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW3M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,2] <- F            # NULL
    ### drop 2
      S <- 0.94 * data[, "FX" ]
      memo <- window( log( S/data[, "FORW3M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW3M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,7] <- F            # NULL
    ### drop 3
      S <- 0.93 * data[, "FX" ]
      memo <- window( log( S/data[, "FORW3M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW3M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,12] <- F           # 28 < 0
    ### rise 1
      S <- 1.05 * data[, "FX" ]
      memo <- window( log( S/data[, "FORW3M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW3M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,17] <- 1 - F
    ### rise 2
      S <- 1.06 * data[, "FX" ] 
      memo <- window( log( S/data[, "FORW3M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW3M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,22] <- 1 - F            # 22 < 0
    ### rise 3
      S <- 1.07 * data[, "FX" ] 
      memo <- window( log( S/data[, "FORW3M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW3M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,27] <- 1 - F            # 81 < 0
  ### 6 months
      tt <- 6 / 12
      F <- data[,"FORW6M"]
      K_A <- F * exp( 0.5 * ( data[,"ATM6m" ]/100 )^2 * tt )
      IV_A <- data[,"ATM6m" ]
      ### 10-Delta options
      delta <- 0.1
      alpha <- -qnorm( delta * exp( log( 1 + data[,"USLIBOR6m" ] / 100 ) * tt ) )
      IV_C <- data[,"BF106m" ] + data[,"ATM6m" ] + 0.5 * data[,"RR106m" ]
      IV_P <- data[,"BF106m" ] + data[,"ATM6m" ] - 0.5 * data[,"RR106m" ]
      K_10C <- F * exp( 0.5 * ( IV_C/100 )^2 * tt + alpha *  IV_C/100 * sqrt( tt ) )
      K_10P <- F * exp( 0.5 * ( IV_P/100 )^2 * tt - alpha *  IV_P/100 * sqrt( tt ) )
      IV_10C <- IV_C
      IV_10P <- IV_P   
      ### 25-Delta options
      delta <- 0.25
      alpha <- -qnorm( delta * exp( log( 1 + data[,"USLIBOR6m" ] / 100 ) * tt ) )
      IV_C <- data[,"BF256m" ] + data[,"ATM6m" ] + 0.5 * data[,"RR256m" ]
      IV_P <- data[,"BF256m" ] + data[,"ATM6m" ] - 0.5 * data[,"RR256m" ]
      K_25C <- F * exp( 0.5 * ( IV_C/100 )^2 * tt + alpha *  IV_C/100 * sqrt( tt ) )
      K_25P <- F * exp( 0.5 * ( IV_P/100 )^2 * tt - alpha *  IV_P/100 * sqrt( tt ) )  
      IV_25C <- IV_C
      IV_25P <- IV_P  
      ### gsi
      sig <- 0.2 * ( IV_A + IV_10C + IV_10P + IV_25C + IV_25P ) 
      gsi_A <- log( K_A/F ) / ( sig/100 * sqrt( tt ) )
      gsi_10C <- log( K_10C/F ) / ( sig/100 * sqrt( tt ) )    
      gsi_10P <- log( K_10P/F ) / ( sig/100 * sqrt( tt ) ) 
      gsi_25C <- log( K_25C/F ) / ( sig/100 * sqrt( tt ) )    
      gsi_25P <- log( K_25P/F ) / ( sig/100 * sqrt( tt ) ) 
      ### matrix with the data
      M <- matrix( NA, ncol = 11, nrow = length( IV_A ) )
      M[,1] <- IV_A/100
      M[,2] <- gsi_A
      M[,3] <- IV_10C/100
      M[,4] <- gsi_10C
      M[,5] <- IV_10P/100
      M[,6] <- gsi_10P
      M[,7] <- IV_25C/100
      M[,8] <- gsi_25C
      M[,9] <- IV_25P/100
      M[,10] <- gsi_25P
      M[,11] <- as.Date( index( IV_A ) )
      colnames( M ) <- c( "IV_A", "gsi_A", "IV_10C", "gsi_10C", "IV_10P", "gsi_10P", "IV_25C", "gsi_25C", "IV_25P", "gsi_25P", "dates" )
      M <- na.omit(M)
      Par <- matrix( NA, nrow = dim( M )[1], ncol = 3 )
      Par[,1] <- M[,1]
      for( i in 1:dim(M)[1] ){
        Par[i,2:3] <- stats::optim( par = c( 0.1, 0.05 ),
                                    fn = meansq_price_error, data = M[i,], gamma0 = M[i,1],
                                    method = "L-BFGS-B" )$par #,
        #control = list(trace = 5) )
      }
    ### drop 1
      S <- 0.92 * data[, "FX" ] 
      sig_bar <- window( sig, as.Date( M[,11] ) )/100
      memo <- window( log( S/data[, "FORW6M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW6M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,3] <- F            # NULL
    ### drop 2
      S <- 0.91 * data[, "FX" ]
      memo <- window( log( S/data[, "FORW6M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW6M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,8] <- F            # 2 < 0
    ### drop 3
      S <- 0.9 * data[, "FX" ]
      memo <- window( log( S/data[, "FORW6M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW6M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,13] <- F           # 36 < 0
    ### rise 1
      S <- 1.08 * data[, "FX" ]
      memo <- window( log( S/data[, "FORW6M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW6M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,18] <- 1 - F
    ### rise 2
      S <- 1.09 * data[, "FX" ] 
      memo <- window( log( S/data[, "FORW6M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW6M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,23] <- 1 - F       # 34 < 0
    ### rise 3
      S <- 1.1 * data[, "FX" ] 
      memo <- window( log( S/data[, "FORW6M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW6M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,28] <- 1 - F       # 117 < 0
  ### 9 months
      tt <- 9 / 12
      F <- data[,"FORW9M"]
      K_A <- F * exp( 0.5 * ( data[,"ATM9m" ]/100 )^2 * tt )
      IV_A <- data[,"ATM9m" ]
    ### 10-Delta options
      delta <- 0.1
      alpha <- -qnorm( delta * exp( log( 1 + data[,"USLIBOR9m" ] / 100 ) * tt ) )
      IV_C <- data[,"BF109m" ] + data[,"ATM9m" ] + 0.5 * data[,"RR109m" ]
      IV_P <- data[,"BF109m" ] + data[,"ATM9m" ] - 0.5 * data[,"RR109m" ]
      K_10C <- F * exp( 0.5 * ( IV_C/100 )^2 * tt + alpha *  IV_C/100 * sqrt( tt ) )
      K_10P <- F * exp( 0.5 * ( IV_P/100 )^2 * tt - alpha *  IV_P/100 * sqrt( tt ) )
      IV_10C <- IV_C
      IV_10P <- IV_P   
    ### 25-Delta options
      delta <- 0.25
      alpha <- -qnorm( delta * exp( log( 1 + data[,"USLIBOR9m" ] / 100 ) * tt ) )
      IV_C <- data[,"BF259m" ] + data[,"ATM9m" ] + 0.5 * data[,"RR259m" ]
      IV_P <- data[,"BF259m" ] + data[,"ATM9m" ] - 0.5 * data[,"RR259m" ]
      K_25C <- F * exp( 0.5 * ( IV_C/100 )^2 * tt + alpha *  IV_C/100 * sqrt( tt ) )
      K_25P <- F * exp( 0.5 * ( IV_P/100 )^2 * tt - alpha *  IV_P/100 * sqrt( tt ) )  
      IV_25C <- IV_C
      IV_25P <- IV_P  
    ### gsi
      sig <- 0.2 * ( IV_A + IV_10C + IV_10P + IV_25C + IV_25P ) 
      gsi_A <- log( K_A/F ) / ( sig/100 * sqrt( tt ) )
      gsi_10C <- log( K_10C/F ) / ( sig/100 * sqrt( tt ) )    
      gsi_10P <- log( K_10P/F ) / ( sig/100 * sqrt( tt ) ) 
      gsi_25C <- log( K_25C/F ) / ( sig/100 * sqrt( tt ) )    
      gsi_25P <- log( K_25P/F ) / ( sig/100 * sqrt( tt ) ) 
    ### matrix with the data
      M <- matrix( NA, ncol = 11, nrow = length( IV_A ) )
      M[,1] <- IV_A/100
      M[,2] <- gsi_A
      M[,3] <- IV_10C/100
      M[,4] <- gsi_10C
      M[,5] <- IV_10P/100
      M[,6] <- gsi_10P
      M[,7] <- IV_25C/100
      M[,8] <- gsi_25C
      M[,9] <- IV_25P/100
      M[,10] <- gsi_25P
      M[,11] <- as.Date( index( IV_A ) )
      colnames( M ) <- c( "IV_A", "gsi_A", "IV_10C", "gsi_10C", "IV_10P", "gsi_10P", "IV_25C", "gsi_25C", "IV_25P", "gsi_25P", "dates" )
      M <- na.omit(M)
      Par <- matrix( NA, nrow = dim( M )[1], ncol = 3 )
      Par[,1] <- M[,1]
      for( i in 1:dim(M)[1] ){
        Par[i,2:3] <- stats::optim( par = c( 0.1, 0.05 ),
                                    fn = meansq_price_error, data = M[i,], gamma0 = M[i,1],
                                    method = "L-BFGS-B" )$par #,
        #control = list(trace = 5) )
      }
    ### drop 1
      S <- 0.9 * data[, "FX" ] 
      sig_bar <- window( sig, as.Date( M[,11] ) )/100
      memo <- window( log( S/data[, "FORW9M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW9M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,4] <- F 
    ### drop 2
      S <- 0.89 * data[, "FX" ]
      memo <- window( log( S/data[, "FORW9M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW9M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,9] <- F    # 2 <0
    ### drop 3
      S <- 0.88 * data[, "FX" ]
      memo <- window( log( S/data[, "FORW9M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW9M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,14] <- F           # 14 < 0
    ### rise 1
      S <- 1.1 * data[, "FX" ]
      memo <- window( log( S/data[, "FORW9M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW9M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,19] <- 1 - F
    ### rise 2
      S <- 1.11 * data[, "FX" ] 
      memo <- window( log( S/data[, "FORW9M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW9M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,24] <- 1 - F
    ### rise 3
      S <- 1.12 * data[, "FX" ] 
      memo <- window( log( S/data[, "FORW9M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW9M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,29] <- 1 - F      # 60 < 0
  ### 12 months
      tt <- 12 / 12
      F <- data[,"FORW12M"]
      K_A <- F * exp( 0.5 * ( data[,"ATM1y" ]/100 )^2 * tt )
      IV_A <- data[,"ATM1y" ]
      ### 10-Delta options
      delta <- 0.1
      alpha <- -qnorm( delta * exp( log( 1 + data[,"USLIBOR1y" ] / 100 ) * tt ) )
      IV_C <- data[,"BF101y" ] + data[,"ATM1y" ] + 0.5 * data[,"RR101y" ]
      IV_P <- data[,"BF101y" ] + data[,"ATM1y" ] - 0.5 * data[,"RR101y" ]
      K_10C <- F * exp( 0.5 * ( IV_C/100 )^2 * tt + alpha *  IV_C/100 * sqrt( tt ) )
      K_10P <- F * exp( 0.5 * ( IV_P/100 )^2 * tt - alpha *  IV_P/100 * sqrt( tt ) )
      IV_10C <- IV_C
      IV_10P <- IV_P   
      ### 25-Delta options
      delta <- 0.25
      alpha <- -qnorm( delta * exp( log( 1 + data[,"USLIBOR1y" ] / 100 ) * tt ) )
      IV_C <- data[,"BF251y" ] + data[,"ATM1y" ] + 0.5 * data[,"RR251y" ]
      IV_P <- data[,"BF251y" ] + data[,"ATM1y" ] - 0.5 * data[,"RR251y" ]
      K_25C <- F * exp( 0.5 * ( IV_C/100 )^2 * tt + alpha *  IV_C/100 * sqrt( tt ) )
      K_25P <- F * exp( 0.5 * ( IV_P/100 )^2 * tt - alpha *  IV_P/100 * sqrt( tt ) )  
      IV_25C <- IV_C
      IV_25P <- IV_P  
      ### gsi
      sig <- 0.2 * ( IV_A + IV_10C + IV_10P + IV_25C + IV_25P ) 
      gsi_A <- log( K_A/F ) / ( sig/100 * sqrt( tt ) )
      gsi_10C <- log( K_10C/F ) / ( sig/100 * sqrt( tt ) )    
      gsi_10P <- log( K_10P/F ) / ( sig/100 * sqrt( tt ) ) 
      gsi_25C <- log( K_25C/F ) / ( sig/100 * sqrt( tt ) )    
      gsi_25P <- log( K_25P/F ) / ( sig/100 * sqrt( tt ) ) 
      ### matrix with the data
      M <- matrix( NA, ncol = 11, nrow = length( IV_A ) )
      M[,1] <- IV_A/100
      M[,2] <- gsi_A
      M[,3] <- IV_10C/100
      M[,4] <- gsi_10C
      M[,5] <- IV_10P/100
      M[,6] <- gsi_10P
      M[,7] <- IV_25C/100
      M[,8] <- gsi_25C
      M[,9] <- IV_25P/100
      M[,10] <- gsi_25P
      M[,11] <- as.Date( index( IV_A ) )
      colnames( M ) <- c( "IV_A", "gsi_A", "IV_10C", "gsi_10C", "IV_10P", "gsi_10P", "IV_25C", "gsi_25C", "IV_25P", "gsi_25P", "dates" )
      M <- na.omit(M)
      Par <- matrix( NA, nrow = dim( M )[1], ncol = 3 )
      Par[,1] <- M[,1]
      for( i in 1:dim(M)[1] ){
        Par[i,2:3] <- stats::optim( par = c( 0.1, 0.05 ),
                                    fn = meansq_price_error, data = M[i,], gamma0 = M[i,1],
                                    method = "L-BFGS-B" )$par #,
        #control = list(trace = 5) )
      }
    ### drop 1
      S <- 0.9 * data[, "FX" ] 
      sig_bar <- window( sig, as.Date( M[,11] ) )/100
      memo <- window( log( S/data[, "FORW12M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW12M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,5] <- F 
    ### drop 2
      S <- 0.88 * data[, "FX" ]
      memo <- window( log( S/data[, "FORW12M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW12M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,10] <- F       # 3 < 0
    ### drop 3
      S <- 0.86 * data[, "FX" ]
      memo <- window( log( S/data[, "FORW12M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW12M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,15] <- F       # 10 < 0
    ### rise 1
      S <- 1.1 * data[, "FX" ]
      memo <- window( log( S/data[, "FORW12M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW12M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,20] <- 1 - F
    ### rise 2
      S <- 1.12 * data[, "FX" ] 
      memo <- window( log( S/data[, "FORW12M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW12M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,25] <- 1 - F
    ### rise 3
      S <- 1.14 * data[, "FX" ] 
      memo <- window( log( S/data[, "FORW12M"] ), as.Date( M[,11] ) ) / ( sig_bar * sqrt(tt) )
      d <- ( window( log( data[, "FORW12M" ]/S ), as.Date( M[,11] ) ) - IV( Par, memo )^2 * tt/2 ) / ( IV( Par, memo ) * sqrt(tt) )
      F <- ( pnorm( -d ) + dnorm( d ) * Par[,1] / sig_bar * ( Par[,2] + 2 * Par[,3] * memo ) )
      A[,30] <- 1- F # test = # 142 von 4039
  ### save the results  
   A_raw <- A
   A <- ifelse( A_raw < 0, 0, A_raw )
   write.xlsx( data.frame( A ), file = "tail prob (Zhang, Xiong) v4.xlsx" )
   write.xlsx( index( F ), file = "memo_Zhang, Xiong v4.xlsx" )     
   write.xlsx( data.frame( A ), file = "tail prob (Zhang, Xiong) v4.xlsx", rowNames = TRUE, overwrite = TRUE )
   write.xlsx( index( F ), file = "memo_Zhang, Xiong v4.xlsx", rowNames = TRUE, overwrite = TRUE )        
  ### plot
   plot(xts(A[,c( 1, 6, 11, 16, 21, 26 )], order.by = index(F)), main ="One month")
   plot(xts(A[,c( 2, 7, 12, 17, 22, 27 )], order.by = index(F)), main ="Three months")   
   plot(xts(A[,c( 3, 8, 13, 18, 23, 28 )], order.by = index(F)), main ="Six months") 
   plot(xts(A[,c( 4, 9, 14, 19, 24, 29 )], order.by = index(F)), main ="Nine months") 
   plot(xts(A[,c( 5, 10, 15, 20, 25, 30 )], order.by = index(F)), main ="Twelve months") 