### this here is the code of the WP Hertrich, Zimmermann or Chapter 3 of my PhD thesis
###
### 
### to check: graphics 2.tex, graphics 3.tex
### 
###
### change some basic options in R
	options( digits = 3, scipen = 8 )
	options( prompt ="R> ", width = 64)
	options( show.signif.stars = FALSE )
### load necessary libraries
	library( "animation" )
	library( fOptions )
	library( "gdata" )
	library( "hydroTSM" )
	library( "lmtest" )
	library( "NMOF" )
	library( "pracma" )
	library( "sandwich" )
	library( xtable )
	library( "xts" )
### parameters
	linestoskip <- 0
### load the required data files
	### EURCHF TS Impl Vola for D = 25\% Calls with different maturities
		### read the implied volatility TS: Call
			data.file.dir <- "/Users/hrm/Documents/PhD/data/FX Options/EURCHF_25/dynivolC.xls"
			data <- read.xls( data.file.dir, method = "csv", sheet = 1, header = TRUE, skip = linestoskip, stringsAsFactors = FALSE, fill = TRUE, na.strings = "" )
			data[, 1] <- as.Date( data[ , 1], "%d/%m/%Y" )
			for( i in 2:7 ){
				data[ , i ] <- zoo( as.numeric( data[ , i ] ), order.by = data[, 1] )
			}
			dynivolC <- xts( data[,-c(1,8,9)], order.by = data[,"date"] )
			dynivolC <- dynivolC[, -7]
	### EURCHF TS Impl Vola for D = -25\% Puts with different maturities
		### read the implied volatility TS
			data.file.dir <- "/Users/hrm/Documents/PhD/data/FX Options/EURCHF_25/dyn impl vol.xls"
			data <- read.xls( data.file.dir, method = "csv", sheet = 1, header = TRUE, skip = linestoskip, stringsAsFactors = FALSE, fill = TRUE, na.strings = "" )
			data[, 1] <- as.Date( data[ , 1], "%d/%m/%Y" )
			for( i in 2:7 ){
				data[ , i ] <- zoo( as.numeric( data[ , i ] ), order.by = data[, 1] )
			}
			dynivol <- xts( data[,-c(1,8,9)], order.by = data[,"date"] )
			dynivol <- dynivol[,1:6]
	### EURCHF TS Impl Vola for at-the-money option contracts with different maturities
		### read the implied volatility TS
			data.file.dir <- "/Users/hrm/Documents/PhD/data/FX Options/atm/atm2"
			data <- read.xls( data.file.dir, method = "csv", sheet = 1, header = TRUE, skip = linestoskip, stringsAsFactors = FALSE, fill = TRUE, na.strings = "" )
			data[, 1] <- as.Date( data[ , 1], "%d/%m/%Y" )
			for( i in 2:7 ){
				data[ , i ] <- zoo( as.numeric( data[ , i ] ), order.by = data[, 1] )
			}
			dynivolA <- xts( data[,-c(1,8,9)], order.by = data[,"date"] )
			dynivolA <- dynivolA[,1:6]			
	### EUR and CHF Libor
		linestoskip <- 1
		### read the CHF 3mLibor
			data.file.dir <- "/Users/hrm/Documents/PhD/data/Libor/LiborCHF3m.xls"
			data <- read.xls( data.file.dir, method = "csv", sheet = 1, header = TRUE, skip = linestoskip, stringsAsFactors = FALSE, fill = TRUE, na.strings = "" )
			data[, 1] <- as.Date( data[ , 1], "%d/%m/%Y" )
			data[ , "libCHF3m" ] <- zoo( as.numeric( data[ , "libCHF3m" ] ), order.by = data[, 1] )
			libCHF3m <- xts( data[ , "libCHF3m" ], order.by = data[,"date"] )
		### read the EUR 3mLibor
			data.file.dir <- "/Users/hrm/Documents/PhD/data/Libor/LiborEUR3m.xls"
			data <- read.xls( data.file.dir, method = "csv", sheet = 1, header = TRUE, skip = linestoskip, stringsAsFactors = FALSE, fill = TRUE, na.strings = "" )
			data[, 1] <- as.Date( data[ , 1], "%d/%m/%Y" )
			data[ , "libEUR3m" ] <- zoo( as.numeric( data[ , "libEUR3m" ] ), order.by = data[, 1] )
			libEUR3m <- xts( data[ , "libEUR3m" ], order.by = data[,"date"] )
		linestoskip <- 0	
		### read the CHF Libor term structure
			data.file.dir <- "/Users/hrm/Documents/PhD/data/Libor/term structure.xls"
			data <- read.xls( data.file.dir, method = "csv", sheet = 1, header = TRUE, skip = linestoskip, stringsAsFactors = FALSE, fill = TRUE, na.strings = "" )
			data <- data[, 1:10]
			for( i in 1:5 ){
				data[, (2 * i) - 1 ] <- as.Date( data[ , (2 * i) - 1 ], "%d/%m/%Y" )
			}
			SF0001M <- xts( na.omit( zoo( as.numeric( data[ , 2 ] ), order.by = data[, 1 ] ) ) )
			SF0002M <- xts( na.omit( zoo( as.numeric( data[ , 4 ] ), order.by = data[, 3 ] ) ) )
			SF0003M <- xts( na.omit( zoo( as.numeric( data[ , 6 ] ), order.by = data[, 5 ] ) ) )
			SF0006M <- xts( na.omit( zoo( as.numeric( data[ , 8 ] ), order.by = data[, 7 ] ) ) )
			SF00012M <- xts( na.omit( zoo( as.numeric( data[ , 10 ] ), order.by = data[, 9 ] ) ) )
		### read the EUR Libor term structure
			data <- read.xls( data.file.dir, method = "csv", sheet = 2, header = TRUE, skip = linestoskip, stringsAsFactors = FALSE, fill = TRUE, na.strings = "" )
			data <- data[, 1:12]
			data[, 1 ] <- as.Date( data[ , 1 ], "%d/%m/%Y" )
			for( i in 2:12 ){
				data[, i ] <- xts( na.omit( zoo( as.numeric( data[ , i ] ), order.by = data[, 1 ] ) ) )
			}
			EE0001M <- xts( na.omit( zoo( as.numeric( data[ , "EE0001M" ] ), order.by = data[, 1 ] ) ) )
			EE0002M <- xts( na.omit( zoo( as.numeric( data[ , "EE0002M" ] ), order.by = data[, 1 ] ) ) )
			EE0003M <- xts( na.omit( zoo( as.numeric( data[ , "EE0003M" ] ), order.by = data[, 1 ] ) ) )
			EE0006M <- xts( na.omit( zoo( as.numeric( data[ , "EE0006M" ] ), order.by = data[, 1 ] ) ) )
			EE0012M <- xts( na.omit( zoo( as.numeric( data[ , "EE0012M" ] ), order.by = data[, 1 ] ) ) )
		linestoskip <- 1
	### EURCHF spot rate
		### read the daily spot file
			data.file.dir <- "/Users/hrm/Documents/PhD/data/FX/EURCHFd.xls"
			data <- read.xls( data.file.dir, method = "csv", sheet = 1, header = TRUE, skip = linestoskip, stringsAsFactors = FALSE, fill = TRUE, na.strings = "" )
			colnames( data ) <- list( "date", "spotd" )
			data[, 1] <- as.Date( data[ , 1], "%d/%m/%Y" )
			data[ , 2 ] <- zoo( as.numeric( data[ , 2 ] ), order.by = data[, 1] )
			spotd <- xts( data[, 2] )
			spotd <- na.omit( spotd )
	### EURCHF FX options vola smile: D = -10%, -25% and due to put-call parity for Calls
		linestoskip<-0
### change margin size
	par(mar= c(4, 4.3, 2, 4.3))
#########
#########
######### 
#########
######### PutVb.pdf
	strike <- 1.2
	b <- c( 0.000001, seq( 0.01, to = strike, by = 0.01 ) )
	spot <- 1.21
	dom <- 0.01
	forn <- 0.03
	vola <- 0.15
	tt <- 0.25
	veestr <- function( b ) {   ## Veestraeten (2008)
		# returns function value at b
    		z1 <- ( log ( spot / strike ) + ( dom - forn + vola^2 / 2 ) * tt ) / ( vola * sqrt( tt ) )
			d2 <- z1 - vola * sqrt( tt )
			z2 <- ( log( b^2 / ( spot * strike ) ) + ( dom - forn + vola^2 / 2 ) * tt ) / ( vola * sqrt( tt ) )
			z3 <- ( log( spot / b ) + ( dom - forn + vola^2 / 2 ) * tt ) / ( vola * sqrt( tt ) )
			z4 <- ( log( b / spot ) + ( dom - forn + vola^2 / 2 ) * tt ) / ( vola * sqrt( tt ) )
			theta <- 2 * ( dom - forn  ) / vola^2
			fb <- strike * exp( - dom * tt ) * pnorm( -d2  ) - spot * exp( - forn * tt ) * ( pnorm( z3 ) - pnorm( z1 ) ) - b * exp( - dom * tt ) * pnorm( -z3 + vola * sqrt( tt ) )  + theta^(-1) * 				( b * exp( - dom * tt ) * pnorm( -z3 + vola * sqrt( tt ) ) - spot * exp( - forn * tt ) *(  b / spot )^(1 + theta) * ( pnorm( z4 ) - pnorm( z2 ) ) - strike * exp( - dom * tt ) * 				( b / strike )^(1 - theta) * pnorm( z2 - theta * vola * sqrt( tt ) ) )
			return( c( fb ) )
		} 
	vepr <- veestr( b )
	z1 <- ( log ( spot / strike ) + ( dom - forn + vola^2 / 2 ) * tt ) / ( vola * sqrt( tt ) )
	d2 <- z1 - vola * sqrt( tt )
	bs <- strike * exp( - dom * tt ) * pnorm( -d2  ) - spot * exp( - forn * tt ) * pnorm( -z1 ) 	
############
############
############ Implied floor (sigma = sigmaVV)
############
############
	delta <- -0.25
	### 
		theta <- 2 * ( data[ , "dom"] - data[ , "for"]  ) / data[ , "vol"]^2
		vepr <- strike * exp( - data[, "dom"] * tt ) * pnorm( -d2  ) - data[, "spot"] * exp( - data[, "for"] * tt ) * pnorm( -z1 ) - b * exp( - data[, "dom"] * tt ) *
				pnorm( -z3 + data[ , "vol"] * sqrt( tt ) ) + data[, "spot"] * exp( - data[, "for"] * tt ) * pnorm( -z3 ) + theta^(-1) * ( b * exp( - data[, "dom"] * tt ) *
				pnorm( -z3 + data[ , "vol"] * sqrt( tt ) ) - data[, "spot"] * exp( - data[, "for"] * tt ) *(  b / data[, "spot"] )^(1 + theta) * ( pnorm( z4 ) - pnorm( z2 ) )
				- strike * exp( - data[, "dom"] * tt ) * ( b / strike )^(1 - theta) * pnorm( z2 - theta * data[ , "vol"] * sqrt( tt ) ) )
		veestr1 <- function( b ) {   ## Veestraeten (2008)
		# returns function value and its derivative at b
    			z1 <- ( log( data[ , "spot"][i] / strike[i] ) + ( data[ , "dom"][i] - data[ , "for"][i] + vola^2 / 2 ) * tt ) / ( vola * sqrt( tt ) )
			d2 <- z1 - vola * sqrt( tt )
			z2 <- ( log( b^2 / ( data[ , "spot"][i] * strike[i] ) ) + ( data[ , "dom"][i] - data[ , "for"][i] + vola^2 / 2 ) * tt ) / ( vola * sqrt( tt ) )
			z3 <- ( log( data[ , "spot"][i] / b ) + ( data[ , "dom"][i] - data[ , "for"][i] + vola^2 / 2 ) * tt ) / ( vola * sqrt( tt ) )
			z4 <- ( log( b / data[ , "spot"][i] ) + ( data[ , "dom"][i] - data[ , "for"][i] + vola^2 / 2 ) * tt ) / ( vola * sqrt( tt ) )
			theta <- 2 * ( data[ , "dom"][i] - data[ , "for"][i]  ) / vola^2
			fb <- strike[i] * exp( - data[, "dom"][i] * tt ) * pnorm( -d2  ) - data[, "spot"][i] * exp( - data[, "for"][i] * tt ) * pnorm( -z1 ) - b * exp( - data[, "dom"][i] * tt ) *
				pnorm( -z3 + vola * sqrt( tt ) ) + data[, "spot"][i] * exp( - data[, "for"][i] * tt ) * pnorm( -z3 ) + theta^(-1) * ( b * exp( - data[, "dom"][i] * tt ) *
				pnorm( -z3 + vola * sqrt( tt ) ) - data[, "spot"][i] * exp( - data[, "for"][i] * tt ) *(  b / data[, "spot"][i] )^(1 + theta) * ( pnorm( z4 ) - pnorm( z2 ) )
				- strike[i] * exp( - data[, "dom"][i] * tt ) * ( b / strike[i] )^(1 - theta) * pnorm( z2 - theta * vola * sqrt( tt ) ) ) - as.numeric( bspr_1[i] )
			dfb <- 1 / ( vola * sqrt( tt )  )  * ( data[, "spot"][i] * exp( - data[, "for"][i] * tt ) / b ) * dnorm( z3 ) + ( 1 - theta ) / theta * exp( - data[, "dom"][i] * tt ) * 
				( pnorm( z4 - theta * vola * sqrt( tt ) ) + dnorm( z4 - theta * vola * sqrt( tt ) ) / ( vola * sqrt( tt ) ) ) - (1 / theta ) * ( b / data[, "spot"][i] )^theta *
				exp( - data[, "for"][i] * tt ) * ( ( 1 + theta ) * ( pnorm( z4 ) - pnorm( z2 ) ) + ( dnorm( z4 ) - 2 * dnorm( z2 ) ) / ( vola * sqrt( tt ) ) ) - (1 / theta ) * ( strike[i] / b )^theta * 
				exp( - data[, "dom"][i] * tt ) * ( ( 1 - theta ) * pnorm( z2 - theta * vola * sqrt( tt ) ) + 2 * dnorm( z2 - theta * vola * sqrt( tt ) ) / ( vola * sqrt( tt ) ) )
			return( c( fb, dfb ) )
		} 
		veestr2 <- function( b ) {   ## Veestraeten (2008)
		# returns function value at b
    			z1 <- ( log( data[ , "spot"][i] / strike[i] ) + ( data[ , "dom"][i] - data[ , "for"][i] + vola^2 / 2 ) * tt ) / ( vola * sqrt( tt ) )
			d2 <- z1 - vola * sqrt( tt )
			z2 <- ( log( b^2 / ( data[ , "spot"][i] * strike[i] ) ) + ( data[ , "dom"][i] - data[ , "for"][i] + vola^2 / 2 ) * tt ) / ( vola * sqrt( tt ) )
			z3 <- ( log( data[ , "spot"][i] / b ) + ( data[ , "dom"][i] - data[ , "for"][i] + vola^2 / 2 ) * tt ) / ( vola * sqrt( tt ) )
			z4 <- ( log( b / data[ , "spot"][i] ) + ( data[ , "dom"][i] - data[ , "for"][i] + vola^2 / 2 ) * tt ) / ( vola * sqrt( tt ) )
			theta <- 2 * ( data[ , "dom"][i] - data[ , "for"][i]  ) / vola^2
			fb <- strike[i] * exp( - data[, "dom"][i] * tt ) * pnorm( -d2  ) - data[, "spot"][i] * exp( - data[, "for"][i] * tt ) * pnorm( -z1 ) - b * exp( - data[, "dom"][i] * tt ) *
				pnorm( -z3 + vola * sqrt( tt ) ) + data[, "spot"][i] * exp( - data[, "for"][i] * tt ) * pnorm( -z3 ) + theta^(-1) * ( b * exp( - data[, "dom"][i] * tt ) *
				pnorm( -z3 + vola * sqrt( tt ) ) - data[, "spot"][i] * exp( - data[, "for"][i] * tt ) *(  b / data[, "spot"][i] )^(1 + theta) * ( pnorm( z4 ) - pnorm( z2 ) )
				- strike[i] * exp( - data[, "dom"][i] * tt ) * ( b / strike[i] )^(1 - theta) * pnorm( z2 - theta * vola * sqrt( tt ) ) ) - as.numeric( bspr_1[i] )
			return( c( fb ) )
		} 
#########	
#########	
######### berechnet implied floor mit VV-Methode
#########	
#########
	### Berechnung der Optionspreise gemss Black-Scholes
		tt <- 3 / 12
		data <- cbind( na.omit( spotd ), log( 1 + na.omit( libCHF3m ) / 100 ), log( 1 + na.omit( libEUR3m ) / 100 ), na.omit( dynivol[,"ivol3m"] ) / 100 )
		data <- na.omit( data )
		colnames( data ) <- list( "spot", "dom", "for", "vol" )
		strike <- data[,"spot"] * exp( qnorm( - delta * exp( data[,"for"] * tt ) ) * data[,"vol"] * sqrt( tt ) + ( data[,"dom"] - data[,"for"] + 0.5 * data[,"vol"]^2 ) * tt )
		bs <- GBSOption(TypeFlag = "p", S = as.numeric( data[ , "spot"] ), X = as.numeric( strike ), Time = 3 / 12, r = as.numeric( data[ ,"dom"] ), 
			b = as.numeric( data[ , "dom" ] - data[, "for"] ), sigma = as.numeric( data[ , "vol" ] ) )
		bspr_1 <- xts( bs@price, order.by = index( data ) )
		plot( bspr_1, minor.ticks = FALSE, auto.grid = FALSE, major.ticks = "months", major.format = "%m.%Y", main = "" )
	### Berechnung der Optionspreise gemss Veestraeten: Vola = VV-Methode
		### Lade die Daten
			data <- cbind( na.omit( spotd ), log( 1 + na.omit( libCHF3m ) / 100 ), log( 1 + na.omit( libEUR3m ) / 100 ), na.omit( dynivol[,"ivol3m"] ) / 100, 
				na.omit( dynivolC[,"ivol3m"] ) / 100, na.omit( dynivolA[,"ivol3m"] ) / 100 )
			data <- na.omit( data )
			colnames( data ) <- list( "spot", "dom", "for", "volP", "volC", "volA" )
		### Berechne Strike fr den 1-M Put-, Call- und ATM-Kontrakt
			strikeP <- data[,"spot"] * exp( qnorm( - delta * exp( data[,"for"] * tt ) ) * data[,"volP"] * sqrt( tt ) + ( data[,"dom"] - data[,"for"] + 0.5 * data[,"volP"]^2 ) * tt )
			strikeC <- data[,"spot"] * exp( -qnorm( - delta * exp( data[,"for"] * tt ) ) * data[,"volC"] * sqrt( tt ) + ( data[,"dom"] - data[,"for"] + 0.5 * data[,"volC"]^2 ) * tt )
			strikeA <- data[,"spot"] * exp( ( data[,"dom"] - data[,"for"] + 0.5 * data[,"volA"]^2 ) * tt )
			strike <- strikeP
			K <- strikeP
			K1 <- lag( strikeP, 1 )
			K2 <- lag( strikeA, 1 )
			K3 <- lag( strikeC, 1 )		
			spot <- lag( data[ , "spot"], 1 )
			dom <- lag( data[ , "dom"], 1 )	
			forn <- lag( data[ , "for"], 1 )		
			volA <- lag( data[,"volA"], 1 )	
			volC <- lag( data[,"volC"], 1 )
			volP <- lag( data[,"volP"], 1 )
		### Berechne D1(K) und D2(K): Wir verwenden fr K den strike price K1 des Folgetages!!!!
			D1 <- ( log( K2 / K ) * log( K3 / K ) ) / ( log( K2 / K1 ) * log( K3 / K1 ) ) * volP
					+ ( log( K / K1 ) * log( K3 / K ) ) / ( log( K2 / K1 ) * log( K3 / K2 ) ) * volA
					+ ( log( K / K1 ) * log( K / K2 ) ) / ( log( K3 / K1 ) * log( K3 / K2 ) ) * volC
					- volA
			d1 <- ( log( spot / K ) + ( dom - forn + volA^2 / 2 ) * tt ) / ( volA * sqrt( tt ) )
			d2 <- d1 - volA * sqrt( tt )
			d1K1 <- ( log( spot / K1 ) + ( dom - forn + volA^2 / 2 ) * tt ) / ( volA * sqrt( tt ) )
			d2K1 <- d1K1 - volA * sqrt( tt )
			d1K3 <- ( log( spot / K3 ) + ( dom - forn + volA^2 / 2 ) * tt ) / ( volA * sqrt( tt ) )
			d2K3 <- d1K3 - volA * sqrt( tt )
			D2 <- ( log( K2 / K ) * log( K3 / K ) ) / ( log( K2 / K1 ) * log( K3 / K1 ) ) * d1K1 * d2K1 * ( volP - volA )^2
					+ ( log( K / K1 ) * log( K / K2 ) ) / ( log( K3 / K1 ) * log( K3 / K2 ) ) * d1K3 * d2K3 * ( volC - volA )^2
			vola <- volA + ( ( - volA + sqrt( volA^2 + d1 * d2 * ( 2 * volA * D1 + D2 ) ) ) / ( d1 * d2 ) )
		### Erweitere die Variable "data"
			data <- cbind( lag( data, 1 ), na.omit( vola ) )
			colnames( data ) <- list( "spot", "dom", "for", "volP", "volC", "volA", "vol" )	
	### 
		theta <- 2 * ( data[ , "dom"] - data[ , "for"]  ) / data[ , "vol"]^2
		vepr <- strike * exp( - data[, "dom"] * tt ) * pnorm( -d2  ) - data[, "spot"] * exp( - data[, "for"] * tt ) * pnorm( -z1 ) - b * exp( - data[, "dom"] * tt ) *
				pnorm( -z3 + data[ , "vol"] * sqrt( tt ) ) + data[, "spot"] * exp( - data[, "for"] * tt ) * pnorm( -z3 ) + theta^(-1) * ( b * exp( - data[, "dom"] * tt ) *
				pnorm( -z3 + data[ , "vol"] * sqrt( tt ) ) - data[, "spot"] * exp( - data[, "for"] * tt ) *(  b / data[, "spot"] )^(1 + theta) * ( pnorm( z4 ) - pnorm( z2 ) )
				- strike * exp( - data[, "dom"] * tt ) * ( b / strike )^(1 - theta) * pnorm( z2 - theta * data[ , "vol"] * sqrt( tt ) ) )
	### given implied vola
		strike <- strike["2011-09-06/"]
		data <- data["2011-09-06/"]
		bspr_1 <- bspr_1["2011-09-06/"]
		vola <- vola["2011-09-06/"]
		theta <- theta["2011-09-06/"]
		k <- dim( data )[1]
		A3m <- data.frame( matrix( NA, nrow = k, ncol = 4 ) )
		colnames( A3m ) <- list( "iter", "estim.prec", "impl b", "Veest Price" )
		for( i in c(1:k) ){
			vola <- data[, "vol"][i] 
			fm1 <- bisect( veestr2, as.numeric( strike[i] - 0.08 ), strike[i], maxiter = 100, tol = .Machine$double.eps^0.5)
			A3m[i,1] <- fm1$iter		# Anzahl Iterationen
			A3m[i,2] <- fm1$estim.prec
			A3m[i,3] <- fm1$root		# geschtztes b
			A3m[i,4]	 <- fm1$f.root		# Veestraetenpreis mit b = root
		}
		Abisect <- A3m
		Abisect <- xts( Abisect, order.by = index( strike ) )
		A3m <- Abisect		
#########	
#########	
pdf( "Fig1.pdf", width = 10, height = 10 )
par( family = "Times" )
	par(mar= c(5.7, 4.7, 2, 4.7))
###	probVS.pdf
	plot( ( spotd - Abisect[,3] )["2011-09-06/"], minor.ticks = FALSE, auto.grid = FALSE, major.ticks = "months", major.format = "%m.%y", main = "", col = 1, 
		ylab = expression( ""*S[t]*" - b, 1.2 - b" ), cex.lab =1.5, las = 2, xaxt = "n", cex.axis = 1.5, lty = "dotted", ylim = c( 0.0042, 0.5 )  )
	drawTimeAxis( ( dynivolC - dynivol )["2010-06-01/"][,1], lab.tstep="quarters", lab.fmt="%m-%Y", las=2, tick.tstep="months" , cex = 1.5)
	points( ( spotd["2011-09-06/"] - Abisect[,3] )["2011-09-06/"], type = "l", lwd = 2, col = "black", lty = "dotted" )	
	points( ( 1.20 - Abisect[,3] )["2011-09-06/"], type = "l", col = "grey65", lwd = 2, pch = 20, lty = "dotted" )
	abline( h = 0 )
	par( new = TRUE )
	plot( Abisect[,3], type = "l", col = "grey", pch = 20, xaxt = "n", yaxt = "n", xlab ="", ylab ="", minor.ticks = FALSE, auto.grid = FALSE, 
		major.ticks = "months", major.format = "%m.%y", ylim = c( 1.05, 1.27 ), main = ""  )
	par( new = TRUE )
	plot( spotd["2011-09-06/2015-01-14"], type = "l", col = "black", pch = 20, xaxt = "n", yaxt = "n", xlab ="", ylab ="", minor.ticks = FALSE, auto.grid = FALSE, 
		major.ticks = "months", major.format = "%m.%y", ylim = c( 1.05, 1.27 ), main = ""  )
		abline( v = attributes(Abisect["2012-09-06"])$index, lwd = 2, col = "grey", lty = 2 )
		text("06.09.2012", y = 0.09, x = attributes(data["2012-12-27"])$index, bty = "n", cex =1.5 )
		abline( v = attributes(Abisect["2012-07-26"])$index, lwd = 2, col = "grey", lty = 2 )
		text("26.07.2012", y = 0.02, x = attributes(data["2012-04-05"])$index, bty = "n", cex =1.5 )
	axis( 4, cex.axis = 1.5, las = 2 )
	abline( h = 1.2 )
	legend( "topright", legend = expression( S[t] - b, 1.20 - b, S[t], b ), lty = c(3,3,1,1), 
		col = c( 1, "grey65", "black", "grey65" ), ncol = 1, lwd = 2, cex = 1.5, bty="n" )
	mtext( expression( ""*S[t]*",  b" ), side=4, las = 3, cex=1.5, line= 3.7 )
par( family = "" )	
dev.off()
#########
#########
#########
#########	
#########	
#########	
#########	
######### RN break probability Veestraeten model (sigmaVV) 1m
#########	
#########
veestrprob <- function( b ) {   ## Veestraeten (2008)
		# returns probability of 
    			z1 <- ( log( data[ , "spot"][i] / 1.20 ) + ( data[ , "dom"][i] - data[ , "for"][i] + vola^2 / 2 ) * tt ) / ( vola * sqrt( tt ) )
			z2 <- ( log( b^2 / ( data[ , "spot"][i] * 1.20 ) ) + ( data[ , "dom"][i] - data[ , "for"][i] + vola^2 / 2 ) * tt ) / ( vola * sqrt( tt ) )
			theta <- 2 * ( data[ , "dom"][i] - data[ , "for"][i]  ) / vola^2
			probV <- pnorm( -z1 + vola * sqrt( tt ) ) - ( b / 1.20 )^( 1 - theta ) * pnorm( z2 - theta * vola * sqrt( tt ) )
			probBS <- pnorm( -z1 + vola * sqrt( tt ) )
			return( c( probV, probBS) )
		} 
	######### RN break probability Veestraeten 1m
		tt <- 1/ 12
		data <- cbind( na.omit( spotd ), log( 1 + na.omit( libCHF3m ) / 100 ), log( 1 + na.omit( libEUR3m ) / 100 ), na.omit( dynivol[,"ivol1m"] ) / 100 )
		data <- na.omit( data )
		data <- data["2011-09-06/2015-01-14"]
		colnames( data ) <- list( "spot", "dom", "for", "vol" )
		B <- data.frame( matrix( NA, nrow = k, ncol = 2 ) )
		colnames( B ) <- list( "probV", "probBS" ) 
		for( i in c(1:k) ){
			vola <- data[, "vol"][i] 
			fm2 <- veestrprob( as.numeric( A1m[i,3] ) )
			# setze negative Werte auf Null
			B[i, 1] <- ifelse( fm2[1] < 0, 0, fm2[1] )
			B[i, 2] <- ifelse( fm2[2] < 0, 0, fm2[2] )
		}	
		veestrprob1m <- B		
		veestrprob1m <- xts( B, order.by = index( data["2011-09-06/" ] ) )
	######### RN break probability Veestraeten 3m
		tt <- 3/ 12
		data <- cbind( na.omit( spotd ), log( 1 + na.omit( libCHF3m ) / 100 ), log( 1 + na.omit( libEUR3m ) / 100 ), na.omit( dynivol[,"ivol3m"] ) / 100 )
		data <- na.omit( data )
		data <- data["2011-09-06/2015-01-14"]
		colnames( data ) <- list( "spot", "dom", "for", "vol" )
		B <- data.frame( matrix( NA, nrow = k, ncol = 2 ) )
		colnames( B ) <- list( "probV", "probBS" ) 
		for( i in c(1:k) ){
			vola <- data[, "vol"][i] 
			fm2 <- veestrprob( as.numeric( A3m[i,3] ) )
			# setze negative Werte auf Null
			B[i, 1] <- ifelse( fm2[1] < 0, 0, fm2[1] )
			B[i, 2] <- ifelse( fm2[2] < 0, 0, fm2[2] )
		}			
		veestrprob3m <- B		
		veestrprob3m <- xts( B, order.by = index( data["2011-09-06/" ] ) )
	######### RN break probability Veestraeten 1y
		tt <- 12 / 12
		data <- cbind( na.omit( spotd ), log( 1 + na.omit( libCHF3m ) / 100 ), log( 1 + na.omit( libEUR3m ) / 100 ), na.omit( dynivol[,"ivol1y"] ) / 100 )
		data <- na.omit( data )
		data <- data["2011-09-06/2015-01-14"]
		colnames( data ) <- list( "spot", "dom", "for", "vol" )
		B <- data.frame( matrix( NA, nrow = k, ncol = 2 ) )
		colnames( B ) <- list( "probV", "probBS" ) 
		for( i in c(1:k) ){
			vola <- data[, "vol"][i] 
			fm2 <- veestrprob( as.numeric( A1y[i,3] ) )
			# setze negative Werte auf Null
			B[i, 1] <- ifelse( fm2[1] < 0, 0, fm2[1] )
			B[i, 2] <- ifelse( fm2[2] < 0, 0, fm2[2] )
		}			
		veestrprob1y <- B		
		veestrprob1y <- xts( B, order.by = index( data["2011-09-06/" ] ) )
#########
#########
######### RN break probability Veestraeten
#########	
#########	
pdf( "Fig2.pdf", width = 10, height = 10 )
par( family = "Times" )
	par(mar= c(5.7, 4.7, 2, 4.7))
###	probVS.pdf
	plot( 100 * na.omit(veestrprob3m[ , 1 ] ), minor.ticks = FALSE, auto.grid = FALSE, major.ticks = "months", major.format = "%m.%y", main = "", col = 1, ylim = c(0, 65 ), 
		ylab = expression("RN Break Probability "*hat(pi)*" (%)" ), cex.lab =1.5, las = 2, xaxt = "n", cex.axis = 1.5 )
	drawTimeAxis( ( dynivolC - dynivol )["2010-06-01/"][,1], lab.tstep="quarters", lab.fmt="%m-%Y", las=2, tick.tstep="months" , cex = 1.5)
	points( 100 * na.omit(veestrprob1m[ , 1 ] ), type = "l", lty = "dotted", lwd = 2, col = "grey30" )	
	points( 100 * na.omit(veestrprob3m[ , 1 ] ), type = "l", col = "black", lwd = 2, pch = 20 )
	points( 100 * na.omit(veestrprob1y[ , 1 ] ), type = "l", lty = "dotted", col = "grey65", lwd = 2, pch = 20 )
	#####
	###### Falls zustzlich der Abstand von 1.20 zu b dargestellt werden soll
	#####
		# points( 100 * ( 1.20 - Abisect[,3] )["2011-09-06/"], type = "l", col = "grey65", lwd = 2, pch = 20, lty = "dotted" )
	legend( "topright", legend = expression( paste( hat( pi )["1m"]^"1.20" ), paste( hat( pi )["3m"]^"1.20" ), paste( hat( pi )["1y"]^"1.20" ) ), lty = c(3, 1, 3), 
		col = c( "grey30", 1, "grey65" ), ncol = 1, lwd = 2, cex = 1.5, bty="n" )
	abline( h = 0 )
		# abline( v = attributes(Abisect["2012-09-06"])$index, lwd = 2, col = "grey", lty = 2 )
		# text("06.09.2012", y = 55, x = attributes(data["2012-12-27"])$index, bty = "n", cex =1.5 )
		# abline( v = attributes(Abisect["2012-07-26"])$index, lwd = 2, col = "grey", lty = 2 )
		# text("26.07.2012", y = 15, x = attributes(data["2012-04-05"])$index, bty = "n", cex =1.5 )
	# axis( 4, at = seq( 0, 60, 10 ), labels = c( "0.0", 0.1, 0.2, 0.3, 0.4, 0.5, 0.6 ), cex.axis = 1.5, las = 2 )
	# mtext("Deviation of FX Rate from Implied Floor", side=4, las = 3, cex=1.5, line= 3.7 )
par( family = "" )	
dev.off()