Strathmore University SU+ @ Strathmore University Library Electronic Theses and Dissertations 2018 Valuation of a locational spread option: the case of tomatoes in Nairobi and Mombasa Counties in Kenya Kenneth Gitonga Kimathi Strathmore Institute of Mathematical Sciences (SIMs) Strathmore University Follow this and additional works at https://su-plus.strathmore.edu/handle/11071/5968 Recommended Citation Kimathi, K. G. (2018). Valuation of a locational spread option: the case of tomatoes in Nairobi and Mombasa Counties in Kenya (Thesis). Strathmore University. Retrieved from http://su- plus.strathmore.edu/handle/11071/5968 This Thesis - Open Access is brought to you for free and open access by DSpace @Strathmore University. It has been accepted for inclusion in Electronic Theses and Dissertations by an authorized administrator of DSpace @Strathmore University. For more information, please contact librarian@strathmore.edu Valuation of a Locational Spread Option: The case of Tomatoes in Nairobi and Mombasa counties in Kenya Kimathi, Kenneth Gitonga Submitted in partial fulfillment of the requirements for the Degree of Master of Science in Mathematical Finance at Strathmore University Institute of Mathematical Sciences Strathmore University Nairobi, Kenya May, 2018 This thesis is available for Library use on the understanding that it is copyright material and that no quotation from the thesis may be published without proper acknowledgement. Declaration I declare that this work has not been previously submitted and approved for the award of a degree by this or any other University. To the best of my knowledge and belief, the thesis contains no material previously published or written by another person except where due reference is made in the thesis itself. c©No part of this thesis may be reproduced without the permission of the author and Strathmore University Kimathi Kenneth Gitonga Signature: ................................................................................................. Date: ......................................................................................................... The thesis of Kimathi Kenneth Gitonga was reviewed and approved by the following: Dr. Samuel Chege Maina Lecturer, Strathmore Institute of Mathematical Sciences Strathmore University Ferdinand Othieno Dean, Strathmore Institute of Mathematical Sciences Strathmore University Prof. Ruth Kiraka Dean, School of Graduate Studies Strathmore University i Abstract Locational Spread Options are financial instruments that can be used by traders wishing to purchase but not physically acquire produce; to hedge their risks, and / or to take speculative positions, based on their knowledge of market dynamics. In this study, we analyze historical tomato price data in Nairobi & Mombasa counties in Kenya; and observe that the Ornstein Uhlenbeck process best captures the price dynamics due to the mean reverting characteristics noted in the deseasonalized price data. We then derive pricing equations and estimate the model parameters via the use of Maximum Likelihood Estimation. Finally, we use these parameter es- timations to perform Monte Carlo simulations, using the antithetic variate variance reduction technique to obtain the option price. ii Contents 1 Introduction 1 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Objectives of the study . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3.1 Main Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3.2 Specific Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.4 Significance of the study . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2 Literature Review 4 3 Research Methodology 7 3.1 The Ornstein - Uhlenbeck model . . . . . . . . . . . . . . . . . . . . . . 7 3.2 Monte Carlo simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 4 Presentation of Research Findings 15 4.1 Parameter Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 4.2 Monte Carlo simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 5 Discussion, Conclusion & Recommendations 25 References 27 Appendix 29 iii List of Figures 4.1 Price & Temperature data plots . . . . . . . . . . . . . . . . . . . . . . . 16 4.2 Deterministic function plots . . . . . . . . . . . . . . . . . . . . . . . . . 18 4.3 Residual regression plots . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.4 At the money spread call option convergence plot . . . . . . . . . . . . . 21 4.5 In the money call option convergence plot . . . . . . . . . . . . . . . . . 22 4.6 Out of the money call option convergence plot . . . . . . . . . . . . . . . 23 iv List of Tables 4.1 Parameter Values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.2 At the money spread call option prices . . . . . . . . . . . . . . . . . . . 20 4.3 In the money spread call option prices . . . . . . . . . . . . . . . . . . . 21 4.4 Out of the money spread call option prices . . . . . . . . . . . . . . . . . 22 4.5 Standard Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4.6 Computation of the Greeks: Delta . . . . . . . . . . . . . . . . . . . . . . 24 v Acknowledgements I would like to thank my supervisor Dr. Samuel Chege Maina, for the guidance and support he offered throughout the writing of this thesis. I am grateful to the Strathmore Institute of Mathematical Sciences (SIMS) for their continued support during my time at Strathmore University. I would also like to thank my colleagues in the MSc. Mathematical Finance Class of 2018, for their support, friendship and encouragement. A special thank you to my friends and family for their prayers and best wishes. Finally, I am deeply thankful to the Lord Almighty for his Amazing Grace. vi Dedication This thesis is dedicated to my daughter, Makena; and to my wife, Karimi; to whom I am thankful for her immense love, support and encouragement. vii 1 Introduction 1.1 Background Agriculture plays a crucial role in Kenya’s economy; as it contributes 26% of the country’s GDP and another 27% indirectly through linages with other economic sectors [1]. Fresh agricultural produce in Kenya is traded on both a wholesale and retail basis, with whole- sale markets located throughout major towns; and the retail markets that are located all over estates and villages. This produce is supplied through a supply chain that consists of farmers, brokers, wholesalers, distributors, retail traders and consumers. Tomatoes are extensively grown for local consumption in Kenya; and can grow in different agro- ecological zones, either under irrigation or rain-fed conditions [2]. Major areas of tomato production include Kirinyaga, Meru Central (Mitunguu & Isiolo), Nyeri, Nakuru (Bahati & Kabazi) and Taita Taveta. In Kenya, commodities are mainly traded via the spot market. Future and forward con- tracts on commodities are yet to be developed for majority of agricultural commodities; but are at an exploratory stage - with futures based on tea being at the forefront of the conversation [3]. The commodity market in Kenya presently consists of the Mom- basa Tea Auction and the Nairobi Coffee Exchange as indicated by the Capital Markets Authority [4]. However, the price information does not reach farmers due to opaque mar- keting systems; to remedy this, the Capital Markets Master Plan states the need for more transparent commodities markets. The plan therefore made a recommendation (amongst others) that a commodity derivatives exchange should be developed for the commodities traded in the spot commodities market. A Spread Option is defined as a financial instrument based upon the difference of the prices of 2 underlying assets. A Locational Spread Option, therefore, is a financial instru- ment based upon price differences between the same commodity in a different location. In the commodity markets, spread options are based on the differences between the prices of the same commodity at 2 different locations (location spreads) or between the prices of the same commodity at 2 different points in time (calendar spreads) or between the prices of inputs to, and outputs from, a production process (processing spreads) as well 1 as between the prices of different grades of the same commodity (quality spreads). In the energy markets, spreads are used to quantify the cost of production of refined products from the raw material used to produce them. A Crack Spread is a simultaneous purchase/sale of crude against the purchase/sale of refined petroleum products. A Spark Spread provides insight into the efficiency of a power plant. It can be defined as the cost of converting a specific fuel (e.g. natural gas) into electricity. In the Agricultural markets, a good example of a spread is the Soybean crush spread traded on the Chicago Board of Trade (CBOT) which comprises of futures contracts of Soybean, Soybean Oil and Soybean meal, which gives an indication of the average gross processing margin, and is used by processors to hedge cash positions, or for pure speculation. The application of commodity spread options in the agricultural markets in Kenya will open up the agro-commodities market to investors seeking a return, as well as traders who may wish to purchase but not physically acquire produce. Therefore, there is a case to be made for the development of financial instruments to be used in the Horticultural sector. This thesis aims to contribute to this development via the pricing of a spread option. 1.2 Problem Statement Derivative contracts on commodities are yet to be developed locally for majority of agri- cultural produce. The Capital Markets Authority has recommended the development of a commodity derivatives exchange, which includes the modelling and pricing of financial instruments to be used by market participants; which are yet to be developed for use in the Kenyan context. 1.3 Objectives of the study 1.3.1 Main Objective 1. To design and implement a locational spread option investigate the price differences that exist for tomatoes in Nairobi and Mombasa counties. 2 1.3.2 Specific Objectives 1. To analyze spot price data of tomatoes in Nairobi and Mombasa, to identify stylized facts such as mean reversion, 2. To ascertain which methodology would suit the development of the model best, 3. To implement the chosen spread option methodology, to develop the locational spread option pricing model, via parameter estimation and application of numerical methods, 4. To compute prices of at the money, in the money and out of the money spread options. 1.4 Significance of the study The application of commodity spread options in the agricultural markets in Kenya will open up the agro-commodities market to investors seeking a return, as well as traders who may wish to purchase but not physically acquire produce. This study therefore aims to contribute to this via modelling and computing the price of a locational spread option. The thesis is arranged as follows: In Chapter 2 we review existing literature on closed form and non-closed form solutions for spread options, in Chapter 3 we discuss the spread option pricing model, which is based on the Ornstein - Uhlenbeck process; as well as the numerical method used (Monte Carlo simulation) to obtain the option price; in Chapter 4, we discuss the parameter estimation results, as well as our numerical pricing results. Finally in Chapter 5, we conclude our study. 3 2 Literature Review Options can either be solved via closed form or non-closed form solutions. Closed form solutions are convenient and allow quick computation of option prices; but may do so at the expense of accuracy and robustness. Non-closed form solutions do not suffer from these disadvantages; but such methods are often computationally involving and time con- suming. Bjerksund & Stensland [5] derive a formula for the spread call value, conditional on fol- lowing the Kirk’s approximation [6], which they show to be a feasible but non-optimal exercise strategy. They then perform numerical investigations, comparing their simula- tions on three models (their model, Kirk’s model [6] and Carmona-Durrleman procedure 7]), using the Monte Carlo simulations as a benchmark; and establish that their for- mula is more precise than the Kirk’s approximation, and marginally more accurate than Carmona-Durrleman’s model. Carmona and Durrleman [7] survey the theoretical and computational problems associ- ated with spread option pricing. They present common features of all the spread options by discussing their roles as speculation devices and risk management tools and review the mathematical framework and numerical algorithms used to price and hedge them. They reviewed a wide scope of existing literature relating to spread option pricing, including closed form pricing models developed by Kirk [6], Bachelier, Samuelson, Carmona & Durrleman [7] amongst others. They outline how pricing and hedging dynamics can be implemented in models for both the spot price and forward curve dynamics. In the Bache- lier model, the underlying indexes are modelled by means of lognormal distributions as prescribed by Samuelson, which is motivated by the desire to reproduce the inherent pos- itivity of the indexes. The positivity restriction does not apply to the spreads themselves. Carmona & Durrleman [7] noted that most pricing algorithms do not address the aspect of hedging, of evaluating the Greeks (partial derivatives of the price). They therefore show that hedging strategies can be computed and implemented in an efficient manner. The closed form formula derived in their paper can be used to compute the Greeks. A comparison of the results of Bachelier’s model with those obtained by the closed form 4 formula derived in [7], and the Kirk’s approximation show that Carmona and Durrle- man’s model is superior to the others as it allows for easy computations of the Greeks. They note that the Geometric Brownian Motion models proposed by Samuelson fail to capture Mean Reversion. This feature is included in historical models by assuming that the dynamics of the underlying indexes Si(t) are given by geometric Ornstein-Uhlenbeck processes instead of Geometric Brownian Motions, which is then represented as a closed form solution in [7]. Commonly used numerical methods to price & hedge financial instruments in the ab- sence of explicit formulae in closed form are as follows: PDE Solvers, Trinomial trees, Monte Carlo computations, and Fourier Transform. Monte Carlo methods give good price approximations but do not address the different sensitivities of the prices [7], which are important as they assist in risk management & hedging. Conversely, Trinomial tree methods allow to compute the partial derivatives along with the price. However, their use results in slow computing times, as well as the fact that they are feasible with 2 assets but may not succeed for higher dimensions. Hurd & Zhou [8] introduce a new formula for general spread option pricing based on Fourier analysis of the payoff function and found it to be stable & applicable in a wide variety of asset pricing models. They point out that where the stock price follows a GBM process and the strike price is greater than 0, there exists a gap in regard to explicit pricing formulae; with several approximation methods available; e.g. numerical integra- tion methods & analytical methods. They provide a numerical integration method for computing spread options in more than 2 dimensions using the FFT. Their method is based on square integrable integral methods for the payoff function, and is applicable to a variety of spread option payoffs in any model for which the characteristic function of the joint return process is given analytically. They showed that FFT provides an accurate and efficient implementation of the pricing formula in low dimensions - but for higher dimensional problems, the issue of dimensionality sets in. Cane and Olivares [9] examine spread option pricing under models with jumps driven by Compound Poission Procesess (CPPs) and stochastic volatilities in the form of Cox 5 Ingersoll Ross (CIR) processes. They derive the Characteristic Functions for 2 market models and use Fast Fourier Transforms (FFTs) to accurately compute spread option prices across a variety of strikes and initial price vectors. They noted that Black & Sc- holes [10] fail to capture critical empirical features of financial markets. They extended Bates’ model [11] to 2 market multivariate models with jumps and stochastic volatil- ity and derived the characteristic function under each model. Using FFT and Hurd & Zhou’s [8] implementation, they produced results which closely matched those produced by Monte Carlo methods in a shorter time. The prices produced were sensitive to jump and correlation parameters. Egorova & Jodas [12] discuss numerical analysis and computing of a spread option pric- ing problem described by a two spatial variable PDE. They develop an explicit difference scheme that retains the benefits of the one-dimensional finite difference method. They then compare these results with those of the Numerical Integration method (NIM) the Fast Fourier Transform (FFT) method & Monte Carlo methods; and establish that the Finite Difference method have the least absolute errors, the computational speed is sec- ond only to that of the NIM, and hence they recommend the use of Finite Difference methods for spread option pricing. Dempster & Hong [13] extended the Fast Fourier Transform technique introduced by Carr & Madan [18] to a multifactor setting for pricing of spread options. They compare the Analytic, Monte Carlo & Fast Fourier Transform methods, and establish that the FFT yields an advantage over the Monte Carlo and Partial Difference methods; mainly due to the reduced computational time of the FFT. Qi Ai [14] checked the exactness of the Kirk [6] and the Bjerksund - Stensland [15] closed form formulae. He found that they presented larger absolute errors when negative strike prices were given; and also proved inadequate in pricing trivariate spread options. He concluded that Numerical Methods such as Monte Carlo methods would provide better results. 6 3 Research Methodology This chapter introduces the model that will be used to model the commodity price, and the pricing technique that will be used to obtain the specific option prices. The determi- nation of the modelling approach, as well as the option pricing technique to be employed is determined from the characteristics of the historical commodity price data. Upon an analysis of our data, we determine that the Ornstein Uhlenbeck process is best suited to model our spread option. We also conclude that Numerical methods, specifically Monte Carlo simulation will be applied to obtain the option price. 3.1 The Ornstein - Uhlenbeck model In most market applications, the underlying indexes are modelled by means of lognormal distributions to reproduce the positivity of the indexes. Therefore, a series of papers proposed to use Arithmetic Brownian Motion (ABM) for the dynamics of spreads [7]. In this way, prices of options can be derived by computing Gaussian integrals leading to simple closed form formulae. This has been found to be inaccurate, as the marginal distribution of the underlying indexes are Gaussian, and can therefore be negative with positive probability, which the ABM cannot accommodate. Therefore, one can assume that the dynamics are given by Geometric Brownian Motion (GBM), which is a convenient basis for the formulation of closed form formulae. The GBM however fails to capture one of the main characteristics of commodity price data - that of mean reversion. An analysis of the residuals derived from historical daily tomato prices in Kenya shows that the prices revert to a long-term mean; which means the GBM is unsuitable in our case. We therefore assume that the dynamics of the underlying indexes in our model follow an Ornstein Uhlenbeck process, which incorporates mean reversion. Alexandridis & Zapranis [16] studied temperature time series data from 7 European cities where Weather Derivatives are traded. They observed that the deseasonalized temperature data residuals follow a mean reverting process, and proceed to model them 7 via an Ornstein Uhlenbeck process, via a model developed by Benth and Saltyte-Benth [17], which follows: Z(t) = µ(t) + (t). (1) where Z(t) is the daily temperature, µ(t) represents the mean process and (t) is the residual process. The mean process is further given as µ(t) = S(t) + p∑ i=1 αi(Z(t− i)− S(t− i)). (2) where S(t) is deterministic and αi, i=1,2,...,p are parameters of the AR(p) process. S(t) plays the role of the long-term average of the temperature, towards which the temperature reverts to. S(t) can therefore be described as the seasonal mean function of temperature. Substituting (2) into (1) yields: Z(t)− S(t) = p∑ i=1 αi(Z(t− i)− S(t− i)) + (t). (3) As long as the residual process has the mean zero, we observe that the expected temper- ature Z˜t = E[Z(t)] follows: Z˜(t)− S(t) = p∑ i=1 αi(Z˜(t− i)− S(t− i)) + (t). (4) The seasonal mean function S(t) is assumed to have the form: S(t) = a0 + a1t+ J∑ j=1 b1j cos(2pij(t− b2j)/365). (5) The dynamics of the deseasonalized temperature are therefore assumed to follow an O-U process as follows: dZ˜(t) = −αZ˜(t)dt+ σ(t)dB(t). (6) where the temperature is defined as Z(t) = S(t) + Z˜(t), α ≥ 0 is a positive constant measuring the speed of mean reversion, and B is a Brownian Motion defined on a filtered probability space (Ω,F ,Ft, P ). 8 Upon analyzing our tomato price data, we noted that it exhibited reversion characteristics as well. We went further to analyze temperature data in Kirinyaga, a major tomato growing region in Kenya [2] and noted that the temperature data reverts to the mean as well. Due to this shared property of mean reversion, we adopt the model developed by [16]& [17] in our case. The price dynamics are given by a Gaussian mean-reverting OU process defined as follows: dT (t) = dS(t)− κ(T (t)− S(t))dt+ σdW (t). (7) where T(t) is the daily price of Tomatoes, κ is the speed of mean reversion, S(t) is a deterministic function modelling trend and seasonality, σ is the volatility of price varia- tions, and W(t) is a Brownian motion. We model S(t) using a Truncated Fourier series to obtain the deterministic function. Its representation is as below: S(t) = a0 + a1(t) + a2cos(2pi(t− a3)/365) + a4(1 + sin(2pi(t− a5)/365)sin(2pit/365). (8) The commodities can then be represented as: dT1(t) = dS1(t)− κ1(T1(t)− S1(t))dt+ σdW 1(t), (9) dT2(t) = dS2(t)− κ2(T2(t)− S2(t))dt+ σdW 2(t). (10) where T1 refers to Tomato prices in Nairobi, and T2 refers to Tomato prices in Mombasa. S1 and S2 are deterministic functions, and W 1 and W 2 are two correlated Brownian Motions with E[dW 1(t)dW 2(t)] = ρij. The prices of Tomatoes in Nairobi and Mombasa are correlated, shown in the matrix below: ∑ =  1 0.6060 0.6060 1  . We therefore convert the vector of correlated Brownian motion processes to a vector of uncorrelated processes: dW 1 dW 2  =  a11 a12 a21 a22 +  dW˜ 1 dW˜ 2  . 9 We select the a′ijs suitable to preserve the correlation structures E(dW i) = 0, E[(dW i)2] = dt, E[dW idW j] = ρijdt ∀i 6= j, which imposes the following conditions: a211 + a 2 12 = 1 a221 + a 2 22 = 1 a11a21 + a21a12 = ρ12 we therefore set a11 = 1, a12 = 0, a21 = ρ12, & a22 = √ 1− ρ212, from which we obtain: dW 1(t) = dW˜ 1(t) dW 2(t) = ρ12dW˜ 1(t) + √ 1− ρ212dW˜ 2(t). Taking T˜ (t) = T(t)-S(t) in equations (9) & (10), we represent the process as follows: dT˜1(t) = −κ1(T˜1(t))dt+ σ1dW 1(t), (11) dT˜2(t) = −κ2(T˜1(t))dt+ σ2dW 2(t). (12) where T˜ (t) = T(t)-S(t). T˜ (t) represents the difference between the tomato price and the deterministic component. Via application of Ito’s lemma, the solution of the Stochastic Differential Equation (SDE) above, T˜ (t) is as follows: T˜1(t) = x1e −κ1t + σ1e−κ1t ∫ t 0 eκ1sdW 1(s), (13) T˜2(t) = x2e −κ2t + σ2e−κ2t ∫ t 0 eκ2sdW 2(s). (14) where x1 = T˜1(0) and x2 = T˜2(0). T˜1(t) and T˜2(t) represent the prices of tomatoes for Nairobi and Mombasa at time t respectively. which further leads to the following pricing equations: T˜1(t) = x1e −κ1t + σ1dW 1(t)(1 + κ), (15) T˜2(t) = x2e −κ2t + σ2dW 2(t)(1 + κ). (16) Proof see Appendix 1. 10 3.2 Monte Carlo simulation This section outlines the pricing method that will be used to establish the option price. Closed form solutions for spread option pricing have been developed for different types of spread options; but should be developed on a case by case basis to avoid modelling errors i.e. a closed form solution based on a model whose underlying indexes follow a GBM process cannot be used for underlying indexes whose deseasonalized prices follow a mean reverting process, as is the case in this study. We have therefore opted to obtain the prices via Monte Carlo simulation, which is performed by generating a large number of sample paths to compute the value of the function of the path whose expectation we evaluate, and averaging these values over the sample paths. Its accuracy be improved by increasing the number of simulations, but this leads to an increase in computation time & cost. However, there are methods that can be used to increase efficiency via the use of reduction of variance techniques. Two of the variation techniques in use are the Anti- thetic variates method and the Control variates method. The antithetic method reduces variance by introducing negative dependence between pairs of replications, whereas con- trol variates take random variables with positive correlation and known expected value under consideration. The spread option SDEs to be simulated are as follows: dT1(t) = dS1(t)− κ1(T1(t)− S1(t))dt+ σdW 1(t), (17) dT2(t) = dS2(t)− κ2(T2(t)− S2(t))dt+ σdW 2(t). (18) where T1 refers to Tomato prices in Nairobi, and T2 refers to Tomato prices in Mombasa. S1 and S2 are deterministic functions, and W 1 and W 2 are two Brownian Motions. The above spread option SDEs will be used to obtain the spread option price. In order to simulate equations (17) & (18) we need to discretize them, and we do this using the Euler discretization scheme. The Euler scheme assumes that 0 = t0 ≤ t1 ≤ 11 ... ≤ tn−1, tn = t and posits the following: dT1(t) = T1(ti+1)− T1(ti), dT2(t) = T2(ti+1)− T2(ti), dS1(t) = S1(ti+1)− S1(ti), dS2(t) = S2(ti+1)− S2(ti), dW 1(t) = √ ti+1 − tiZ1i+1, dW 2(t) = √ ti+1 − tiZ2i+1. therefore, we substitute in equation (17) and (18) to obtain: T1(ti+1)− T1(ti) = S1(ti+1)− S1(ti)− κ1(T1(ti)− S1(ti)(ti+1 − ti + σ1 √ ti+1 − tiZ1i+1,(19) T2(ti+1)− T2(ti) = S2(ti+1)− S2(ti)− κ2(T2(ti)− S2(ti)(ti+1 − ti + σ2 √ ti+1 − tiZ2i+1.(20) Assuming a fixed grid spacing ti+1 − ti = h, then ti = ih. Therefore, equation (19) becomes: T1(i+ 1) = T1(i) + S1(i+ 1)− S1(i)− κ1(T1(i)− S1(i)) ∗ h+ σ1 √ hZ1i+1, (21) T2(i+ 1) = T2(i) + S2(i+ 1)− S2(i)− κ2(T2(i)− S2(i)) ∗ h+ σ2 √ hZ2i+1. (22) The expectation function for spread option pricing is therefore given as below. f(t) = e−rTE[(T2(i+ 1)− T1(i+ 1)−K|Ft)]+. The local truncation error (LTE) of the Euler method is defined as the error made in a single step i.e. the difference between the numerical solution after one step, y1 and the exact solution at time t1 = t0 + h, given by: y1 − y0 = hf(t0, y0), y1 = y0 + hf(t0, y0). For the exact solution, we make use of the Taylor expansion of the function y around t0. y(t0 + h) = y(t0) + hy ′(t0) + 1 2 h2y”(t0) +O(h 3). The LTE is given by the difference between these equations: y0(t0 + h)− y1 = 1 2 h2y′′(t0) +O(h3). The algorithm to be used is as follows: 12 1. Fix the number of monitoring dates & time step ∆ = T/N . 2. Starting from T (0) = x0, simulate the spot prices T1(i + 1) and T2(i + 1) as given in equations (21) & (22). The increment W j(t∆) −W j((t − 1)∆) (where j=1,2) is simulated according to a normal distribution N(0,∆): W (j)(t∆)−W (j)((t− 1)∆) = √ ∆ ∗ φ−1(u(j)i ). (23) where u is a uniform (0,1) random variable and j−1 is the inverse cumulative dis- tribution of the standard normal distribution. 3. Update the average according to Aj(i∆) = i− 1 i xAj((i− 1)∆) + S(i∆) i , (24) Aj(0) = s0. (25) where j=1,2. 4. Compute the discounted option payoff: f(t) = e−rTE[(T2(i+ 1)− T1(i+ 1)−K|Ft)]+. 5. Repeat step 2 to 4 & discount to obtain the option price and average across simu- lations: p˜ = 1 N ∑ p. (26) 6. Evaluate accuracy by computing the standard error: se = √ σ2 N . (27) 7. The confidence interval is given by: p˜+ Z1−α/2 ∗ se. (28) The antithetic variates method was used. This was done via the following steps: 13 (a) Computing a vector U of Normal random variables, (b) Computing a vector of past prices S1 under the O-U process, (c) Creating a new vector V given by V=-U, (d) Computing a new vector of asset prices S2 under the O-U process based on vector V, (e) Taking an average of all the values consisting both vectors S1 and S2 that we call A1 and A2, (f) Discounting the average values obtained to obtain the option price. The Greeks measure sensitivity of the prices and are used by practitioners to manage risk. There are 4 major Greeks [20]: • Delta, which measures the rate of change of the option price with respect to the price of the underlying asset, • Gamma, which is the rate of change of the Delta with respect to the price of the underlying asset, • Theta, which measures the rate of change of the option value with respect to the passage of time, and • Vega, which is the rate of change of the value of the portfolio with respect to the volatility of the underlying asset. In order to numerically compute the Delta via Monte Carlo simulations, we proceed to compute the limit as ξ → 0 of the following expression [7]: ∆1(ξ) ≈ 1 ξ (P (x1 + ξ)− P (x1)). ∆1(ξ) is computed for each ξ in a sequence going to zero. The simulation is conducted for x1 and again for x1 + ξ [7]. 14 4 Presentation of Research Findings In this chapter we analyze past tomato price and local temperature data to review data attributes that would inform our model choice. We establish that the price & temper- ature data exhibit mean reverting properties; therefore confirming the appropriateness of the Ornstein Uhlenbeck model we adopted. Additionally, we obtain parameters from our historical price data; and proceed to use them to price the option via Monte Carlo simulations. We conclude by discussing pricing results for call spread options for At the Money (ATM), In the Money (ITM) and Out of the Money (OTM) options. 4.1 Parameter Estimation Our market model is based on the dynamics of daily spot commodity prices, and is there- fore reliant on historical price data to inform our modelling approach and parameter estimation.1 Commodity price data for the period February 2014 to February 2018 for a crate (64kg) of tomatoes in Nairobi and Mombasa was obtained and used for purposes of our study. The price data comprises of weekday market prices, and hence weekends are not included in our study. Data outliers were removed using the 3 sigma method, which proved more effective than the interquartile range method (IQR), which when used did not remove all outliers. Moreover, the data contained some missing values, which were filled in using linear interpolation. Additionally, daily mean temperature readings2 for Kirinyaga county (a major tomato growing region in Kenya) for the period February 2014 to February 2018 have been ob- tained and analyzed. The analysis of this particular data set was done only to inform the choice of model; since temperature and price data share the stylized fact of mean reversion [16]. Therefore, the parameters from this data have not been used in the computation of our option price. The raw price data for tomato prices in Nairobi & Mombasa, as well as that for temperature is shown in Figure 1. A cursory look at the data shows that in 1The data used in this study is derived from local daily commodity prices compiled by the Ministry of Agriculture, published in the Business Daily newspaper on each weekday. 2This data was obtained from the Climate Forecast system, which is availed by the National Centres for Environmental Protection, based in the USA. 15 all cases, the data is characterized by cycles that exhibit annual seasonality cycles; with some peaks and troughs specifically in the price data. 0 200 400 600 800 1000 1200 4000 6000 8000 10000 Pr ic es in K ES Tomato Price data - Nairobi Tomato Price 0 200 400 600 800 1000 1200 0 5000 10000 Pr ic es in K ES Tomato Price data - Mombasa Tomato Price 0 200 400 600 800 1000 1200 10 15 20 25 ° ce n tig ra de Temperature data - Kirinyaga Temperature Figure 4.1: Price data in Nairobi, Mombasa & Temperature data in Kirinyaga Recall the O-U model for the price dynamics given in equation (7), which includes a deterministic component, S(t), which models trend & seasonality. We begin by modelling this component for both the price and temperature data, using the Truncated Fourier series given in equation (8). The graphical representation is shown in Figure 2, with the parameter values of S(t) given in Table 1. We then determine the residuals by differencing the deterministic component from the raw data in each case (i.e. for price and temperature data). A multiple linear regression was then applied to the residuals to obtain the residual regression plots as represented in Figure 3. In all cases, the data is seen to revert to a long-term mean. Having confirmed that the price residuals are mean reverting, we proceed to estimate the speed of mean reversion, κ, as well as the standard deviation, σ. We proceed to obtain these parameters via discretization of the O-U process 16 as follows: ∆T1(t) = ∆S1(t) + κ1(T1(t− 1)− S1(t− 1)) + σ1 √ ∆tZ1(t), (29) ∆T2(t) = ∆S2(t) + κ2(T2(t− 1)− S2(t− 1)) + σ2 √ ∆tZ2(t). (30) Expanding the above equation with ∆t = 1, we have T (t)− T (t− 1) = S(t)− S(t− 1) + κ(T (t− 1)− S(t− 1)) + σ(t) (31) this arises from the fact that ∆T (t) = T (t)− T (t− 1), ∆S(t) = S(t)− S(t− 1), (t) = √ ∆tZ(t) By rearranging, we have that T (t)− S(t) = T (t− 1)− S(t− 1) + κ(T (t− 1)− S(t− 1)) + σ(t). (32) We then set T˜ (t) = T (t)− s(t) to give: T˜ (t) = T˜ (t− 1) + κT˜ (t− 1) + σ(t). (33) Equivalently, T˜ (t) = (1 + κ)T˜ (t− 1) + σ(t). (34) Substituting with α = 1 + κ, our model is reduced to: T˜ (t) = αT˜ (t− 1) + σ(t). (35) Equation (33) is a simple AR(1) model, hence the parameters κ and σ are estimated via the Maximum Likelihood Estimation method (MLE). 17 0 200 400 600 800 1000 1200 4000 6000 8000 10000 Pr ic es in K ES Tomato Price data - Nairobi Tomato Price Deterministic function 0 200 400 600 800 1000 1200 0 5000 10000 Pr ic es in K ES Tomato Price data - Mombasa Tomato Price Deterministic function 0 200 400 600 800 1000 1200 10 15 20 25 ° ce n tig ra de Truncated Fourier Transform - Temperature data Temperature Deterministic function Figure 4.2: Truncated fourier transform - Price & Temperature data. This shows the deterministic (deseasonalized) component, S(t) of our model. 18 0 200 400 600 800 1000 1200 -2000 -1000 0 1000 Pr ic es in K ES Residual Regression - Nairobi data Residuals 0 200 400 600 800 1000 1200 -2000 0 2000 4000 Pr ic es in K ES Residual Regression - Mombasa data Residuals 0 200 400 600 800 1000 1200 -5 0 5 ° ce n tig ra de Residual Regression - Temperature data Residuals Figure 4.3: Residual regression - Price & Temperature data. These residuals were ob- tained from differencing the deterministic component from the raw data. We observe that the price and temperature data reverts to the mean. We obtained the following model parameters for Nairobi & Mombasa: Table 4.1: Parameter Values Parameter Nairobi Mombasa a0 4.9743 5.9010 a1 -0.0001 -0.0007 a2 -0.8866 0.7864 a3 -2.4918 -0.1943 a4 0.7568 0.2999 a5 -0.0273 0.0297 σ 254.1715 401.4495 κ 0.9610 0.9639 The parameters a0, a1, a2, a3, a4&a5 were estimated from the deseasonalized component 19 S(t) represented in equation (8) via the Truncated Fourier Transform. The parameters σ and κ were estimated via Maximum Likelihood Estimation on the residual data, following the discretization of equation (11) & (12) into an AR(1) process represented by equation (35). The values of σ for Nairobi and Mombasa are noted to be quite high at 254.1715 and 401.4495 respectively, and the speed of mean reversion κ for both the residuals for Nairobi and Mombasa price data are very close to 1, showing that the residuals revert to the mean at a high rate. These parameters were used to perform the Monte Carlo simulation described in the following section. 4.2 Monte Carlo simulation The Monte Carlo simulation will enable us to price our spread option. In order to address the efficiency & computing time problem, we used the Antithetic variate variance reduction technique. We simulated call prices for at the money (ATM), in the money (ITM) and out of the money (OTM) options under different number of simulations, results of which are displayed below. We performed these simulations using an Intel (R) Core (TM) I5-6267U CPU, with a clock size of 2.90 GHz. No. of Simulations Option Price Computation time (seconds) 10 0 0.046 100 12.1750 0.034 1,000 9.3936 0.071 10,000 10.1334 0.378 100,000 9.4405 3.789 1,000,000 9.4792 47.140 10,000,000 9.5861 729.720 Table 4.2: ATM Spread call option prices computed with Monte Carlo simulations. Pa- rameter values are S1(t) = 6000, S2(t) = 5000, σ1 = 254.1715, σ2 = 401.4495, Strike Price: K=1000, Time to maturity: T=1, Risk free rate: r=0.07. 20 0 2 4 6 8 10 No. of simulations 106 0 2 4 6 8 10 12 14 O pt io n Pr ice in K ES 0 100 200 300 400 500 600 700 800 Ti m e in s ec on ds ATM spread call option prices Figure 4.4: At the money spread call option prices. We observe that the option prices converge after approximately 100,000 simulations. No. of Simulations Option Price Computation time (seconds) 10 94.8737 0.037 100 75.9642 0.041 1,000 69.7833 0.051 10,000 61.8587 0.677 100,000 62.1036 8.339 1,000,000 62.1976 91.516 10,000,000 62.7888 720.458 Table 4.3: ITM Spread call option prices computed with Monte Carlo simulations. Pa- rameter values are S1(t) = 6000, S2(t) = 5000, σ1 = 254.1715, σ2 = 401.4495, Strike Price: K=500, Time to maturity: T=1, Risk free rate: r=0.07. 21 0 2 4 6 8 10 No. of simulations 106 60 65 70 75 80 85 90 95 O pt io n Pr ice in K ES 0 100 200 300 400 500 600 700 800 Ti m e in s ec on ds ITM spread call option prices Figure 4.5: In the money spread call option prices. We observe that the option prices converge after approximately 100,000 simulations. No. of Simulations Option Price Computation time (seconds) 10 0 0.040 100 0 0.036 1,000 0 0.0047 10,000 0.0472 0.720 100,000 0.0362 8.907 1,000,000 0.0316 101.538 10,000,000 0.0316 1127.878 Table 4.4: OTM Spread call option prices computed with Monte Carlo simulations. Pa- rameter values are: S1(t) = 6000, S2(t) = 5000, σ1 = 254.1715, σ2 = 401.4495, Strike Price: K=2000, Time to maturity: T=1, Risk free rate: r=0.07. 22 0 2 4 6 8 10 No. of simulations 106 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 O pt io n Pr ice in K ES 0 200 400 600 800 1000 1200 Ti m e in s ec on ds OTM spread call option prices Figure 4.6: Out of the money spread call option prices. We observe that the option prices converge after approximately 100,000 simulations. No. of Simulations Standard Error - Nairobi Standard Error - Mombasa 10 80.3761 126.9495 100 25.4172 40.1450 1,000 8.0376 12.6949 10,000 2.5417 4.0145 100,000 0.8038 1.2695 1,000,000 0.2542 0.4014 10,000,000 0.0804 0.1269 Table 4.5: Standard Errors computed with Monte Carlo simulations. In Table 2, 3 & 4 we show the results from a Monte Carlo simulation performed for spread call options at the money, in the money and out of the money. The option prices for in the money options are of course much higher as the writer requires a higher premium. The 23 converse is true for out of the money options. Figures 4, 5 & 6 graphically represent the speed of convergence and computation times. We established that spread option prices in all cases (ATM, ITM and OTM) converge after approximately 100,000 simulations. Additionally, Table 5 shows that the standard error at this convergence point is also seen to be low in comparison to those seen for fewer simulations. We therefore conclude that the appropriate number of simulations to estimate the spread option price is 100,000. We observe that the computation time in this case is relatively low (3.789 sec for ATM spread call options, 8.339 sec for ITM spread call options and 8.907 for OTM spread call options); as compared to 729.720 sec, 720.458 sec and 1127.878 sec for 10,000,000 simulations for ATM, ITM and OTM spread call options respectively) which therefore presents a good balance between accuracy and computation time. Additionally, we simulated the value of the Greeks i.e. Delta via Monte Carlo simulation, as shown below: ξ ∆1(ξ) 10 0.01562 9 0.028177778 8 0.0064625 7 -0.002728571 6 0.0126 5 -0.00946 4 0.0066 3 0.0576 2 0.08005 1 0.1971 Table 4.6: Delta computed with Monte Carlo simulations, by re-computing the option price with a slight change in the underlying price, S1 = 6000. The average delta in this case comes to 0.03920, which means that when the price of tomatoes changes by a small amount, the price of tomatoes changes by about 3.92% of this amount. 24 5 Discussion, Conclusion & Recommendations We set out to value a locational spread option on Tomatoes in Nairobi and Mombasa counties in Kenya; and did so by determining an appropriate model based on the char- acteristics of historical data, estimated the parameters of the model, and proceeded to apply numerical methods to determine the price of the spread option. We have demon- strated that an Ornstein Uhlenbeck process is suitable for modelling the locational spread option; as it incorporates mean reversion characteristics of the historical prices. More- over, we obtained the option price via the use of Monte Carlo simulations, in particular via use of the antithetic variates method. In our approach, we analyzed historical price data to determine whether the GBM, widely used in previous literature, could be used to model the spread in our case. We found that the prices were mean reverting, immediately disqualifying the GBM. We showed that an O-U model was better placed to model the option. Our price dynamics hence follow a Gaussian mean - reverting process as given in equations (9 & 10). We then proceed to model the deterministic function, S(t), which is given in equation (8) and features in equations (9 & 10); and thereafter applied Ito’s lemma to obtain the pricing equations for Tomato prices in Nairobi and Mombasa. We then estimated the parameters via MLE, using the pricing equations given in equa- tion (13 & 14), and subsequently performed Monte Carlo simulations; making use of the antithetic variate technique to hasten the convergence of the prices. We obtained option prices for a varied number of simulations for at the money, in the money and out of the money spread options, and summarized these results in Table 2, 3 & 4. We showed that generally all these options converged after 100,000 Monte Carlo simulations; as repre- sented graphically in figure 4, 5 & 6. We therefore conclude that 100,000 Monte Carlo simulations are appropriate to compute the spread option price, as they provide a good compromise in terms of speed of convergence and pricing accuracy. In future research, the model can be modified to compute the speed of mean reversion not as a constant, but as a function of time (i.e. κ(t) as opposed to κ), as it can be argued that the speed of mean reversion is not constant for the entire time period in the data set, as shown in Alexandridis & Zapranis [16]. Moreover, further research can model volatility as a seasonal function, in order cater for the seasonality observed in the price residuals as 25 observed in Benth & Saltyte-Benth [19]. Future research can be conducted to perform a comparison of numerical pricing methods for this spread option, similar to the work done by Egorova & Jodas [12]; to compare the results from our Monte Carlo simulation to other numerical methods such as Finite Difference, Fast Fourier Transform & Numerical Integration method, to establish which of these methods gives a better computation of the spread option price, from an accuracy, robustness and computational efficiency point of view. More research can be done in modelling the deterministic component of the spread option pricing model; as it can also be modelled using Wavelet Analysis as done by Alexandridis & Zapranis [16], instead of the Truncated Fourier transform used in our discussion. Finally, further research could be conducted on this spread option pricing problem, to develop a closed form solution gives us the price of the option, as well as provides computations for the Greeks, which have proven useful in other markets for hedging purposes. 26 References 1. Food and Agriculture Organization of the United States. Kenya at a glance. Re- trieved from http://www.fao.org/kenya/fao-in-kenya/kenya-at-a-glance/en/ 2. Larsen, K., Kim, R., & Theus, F. (Eds.). (2009). Agribusiness and innovation systems in Africa. World Bank Publications. 3. Adelaide Changole. (2016, June 16). [Worlds First Tea Futures Contracts May Be Started in Kenya.Bloomerg Terminal Retrieved from https://www.bloomberg.com/ news/articles/2016-06-15/world-s-first-tea-futures-contracts-may-be-introduced-in-kenya 4. The Capital Markets Authority (2016). The development of commodity derivatives markets in Africa [PDF slides]. Retrieved from United Nations Conference on Trade and Development (UNCTAD) website: http://unctad14.org/Documents/gcf2016 Contribution Kenya CommodityExchange en.pdf 5. Peter Bjerksund & Gunnar Stensland (2006). Closed Form Spread Option Valua- tion. 6. E. Kirk (1995). Correlation in the Energy Markets, in Managing Energy Price Risk. London: Risk Publications and Enron. 7. Rene Carmona & Valdo Durrleman (2003). Pricing and Hedging Spread Options. 8. Hurd, T. R., & Zhou, Z. (2010). A Fourier transform method for spread option pricing. SIAM Journal on Financial Mathematics, 1(1), 142-157. 9. Olivares, P., & Cane, M. (2014). Pricing Spread Options under Stochastic Corre- lation and Jump-Diffusion Models. arXiv preprint arXiv:1409.1175. 10. Black, F., & Scholes, M. (1973). The pricing of options and corporate liabilities. Journal of political economy, 81(3), 637-654. 11. Bates, D. S. (1996). Jumps and stochastic volatility: Exchange rate processes implicit in deutsche mark options. The Review of Financial Studies, 9(1), 69-107. 12. Egorova, V. N., & Jdar, L. (2016). An efficient method for solving spread op- tion pricing problem: numerical analysis and computing. In Abstract and Applied Analysis (Vol. 2016). Hindawi. 27 13. Dempster, M. A. H., & Hong, S. G. (2002). Spread option valuation and the fast Fourier transform. In Mathematical FinanceBachelier Congress 2000 (pp. 203- 220). Springer, Berlin, Heidelberg. 14. Ai, Q. (2014). Pricing of Spread Options in Energy Markets with Non-zero Strikes (Master’s thesis). 15. Bjerksund, P., Stensland, G., & Vagstad, F. (2011). Gas storage valuation: Price modelling v. optimization methods. The Energy Journal, 203-227. 16. Alexandridis, A., & Zapranis, A. D. (2012). Weather derivatives: modeling and pricing weather-related risk. Springer Science & Business Media. 17. Benth, J. S., & Benth, F. E. (2012). A critical view on temperature modelling for application in weather derivatives markets. Energy Economics, 34(2), 592-602. 18. Carr, P., & Madan, D. (1999). Option valuation using the fast Fourier transform. Journal of computational finance, 2(4), 61-73. 19. Benth, F. E., & Benth, J. S. (2007). The volatility of temperature and pricing of weather derivatives. Quantitative Finance, 7(5), 553-561. 20. Global Association of Risk Professionals. (2015). Financial Risk Manager, FRM Exam Part I, Valuation and Risk Models. In The Greek Letters (pp. 145-164), Pearson. 28 Appendix Appendix 1 Solution of SDE of T˜ via Ito’s lemma dT (t) = dS(t)− κ(T (t)− S(t))dt+ σdW (t). (36) Via Ito’s lemma, dy = [ ∂y ∂t + µ(xt) ∂y ∂x + 1 2 σ2(xt) ∂2y ∂x2 ] dt+ σ(xt) ∂y ∂x dW, dy = [ κT˜ eκt − κT˜ eκt + 1 2 σ(0) ] dt+ σeκtdW,∫ t 0 d(eκsT˜ ) = ∫ t 0 σeκsdW, eκtT˜t − eκtT˜0 = ∫ t 0 σeκsdWs, T˜t = T˜0e −κt + σe−κt ∫ t 0 eκsdWs. Through integration by parts, i.e.∫ udv = uv − ∫ vdu. We evaluate the integral in the above equation as follows:∫ t 0 eκsdWs = e κsWs|t0 −Wsκeκs|t0, = [eκtWt − eκ0W0]− [Wtkeκt −W0κek0], = eκtzt −W0 −Wtκeκt +W0κ, = Wte κt(1 + κ). Therefore T˜t = T˜0e −κt + σe−kt[Wtekt(1 + κ)], (37) T˜t = T˜0e −κt + σ[Wt(1 + κ)]. (38) The price dynamics for Tomato prices in Nairobi and Mombasa are as below. T˜1(t) = T˜1(0)e −κ1t + σ1[W1(t)(1 + κ1)], (39) T˜2(t) = T˜2(0)e −κ2t + σ2[W2(t)(1 + κ2)]. (40) 29 Appendix 2 MATLAB code 1 -----Locational Spread Option Pricing MATLAB Code----- 2 --------------------------Data read-------------------------- 3 4 clear; clc; close all; 5 6 B = xlsread(’Tomato_det_int.xlsx’, ’Sheet1’); 7 8 9 Date = B(:,1); 10 Nairobi = B(:,2); 11 Mombasa = B(:,3); 12 13 MATLABDate = x2mdate(Date,0,’datetime’) 14 15 save(’MATLABDate.mat’) 16 save(’Nairobi.mat’) 17 save(’Mombasa.mat’) 18 --------------------------TruncatedF-------------------------- 19 20 function Gamma = TruncatedF(a0,a1,a2,a3,a4,a5,t) 21 % Seasonal function with trend, the truncated Fourier transform 22 %Gamma = TruncatedF(x(1), x(2), x(3), x(4), x(5), x(6), t); 23 Gamma = a0 + a1*t + a2*cos(2*pi*(t-a3)/365)+ a4*(1+sin(2*pi*(t - a5)/365)).*sin(2*pi*t/365); 24 end 25 26 --------------------------TruncatedFmsa-------------------------- 27 28 function Beta = TruncatedF(a0,a1,a2,a3,a4,a5,t) 29 % Seasonal function with trend, the truncated Fourier transform 30 %Gamma = TruncatedF(x(1), x(2), x(3), x(4), x(5), x(6), t); 31 Beta = a0 + a1*t + a2*cos(2*pi*(t-a3)/365)+ a4*(1+sin(2*pi*(t - a5)/365)).*sin(2*pi*t/365); 32 end 33 34 --------------------------Est_Sim_Msa-------------------------- 35 36 %clear; clc; close all; 37 38 load(’Mombasa.mat’) 39 y = length(Mombasa); 40 t = (1:1:y)’; 41 42 init = [20, 1, 1, 10, 1, 10]; 43 options = optimoptions(@lsqnonlin,’Algorithm’,’trust-region-reflective’); 44 modelfun = @(b) b(1) + b(2)*t + b(3)*cos(2*pi*(t - b(4))/365)... 45 + b(5)*(1 + sin(2*pi*(t - b(6))/365)).*sin(2*pi*t/365) - Mombasa; 46 y = lsqnonlin(modelfun,init,[],[],options) 47 Beta = TruncatedFmsa(y(1), y(2), y(3), y(4), y(5), y(6), t); 48 49 M=Mombasa-Beta 50 51 M 52 53 [B,BINT,W2] = regress(M(2:end),M(1:end-1)) 54 30 55 -----Est_Sim_Nbi----- 56 57 %clear; clc; close all; 58 59 load(’Nairobi.mat’) 60 x = length(Nairobi); 61 t = (1:1:x)’; 62 63 init = [20, 1, 1, 10, 1, 10]; 64 options = optimoptions(@lsqnonlin,’Algorithm’,’trust-region-reflective’); 65 modelfun = @(b) b(1) + b(2)*t + b(3)*cos(2*pi*(t - b(4))/365)... 66 + b(5)*(1 + sin(2*pi*(t - b(6))/365)).*sin(2*pi*t/365) - Nairobi; 67 x = lsqnonlin(modelfun,init,[],[],options) 68 Gamma = TruncatedF(x(1), x(2), x(3), x(4), x(5), x(6), t); 69 70 J=Nairobi-Gamma 71 72 J 73 74 [B,BINT,W1] = regress(J(2:end),J(1:end-1)) 75 76 77 --------------------------KappaNairobi-------------------------- 78 % Prices at t, X(t) 79 Pt = J(2:end); 80 81 % Prices at t-1, X(t-1) 82 Pt_1 = J(1:end-1); 83 84 % Discretization for daily prices 85 dt = 1/365; 86 87 % PDF for discretized model 88 mrjpdf = @(Pt, phi, sigmaSq) ... 89 exp((-(Pt-phi.*Pt_1).^2)/(2.*sigmaSq)).* ... 90 (1/sqrt(2.*pi.*sigmaSq)); 91 92 % Constraints: 93 lb = [-Inf 0]; 94 ub = [1 Inf]; 95 96 % Initial values 97 x0 = [0 var(J)]; 98 99 % Maximum likelihood estimation 100 params_Nbi = mle(Pt,’pdf’,mrjpdf,’start’,x0,’lowerbound’,lb,’upperbound’,ub,... 101 ’optimfun’,’fmincon’); 102 103 % Calibrated parameters 104 105 kappa_Nai = params_Nbi(1) 106 107 sigma_Nai = sqrt(params_Nbi(2)) 108 109 -----KappaMombasa----- 110 111 % Prices at t, X(t) 112 Pt = M(2:end); 31 113 114 % Prices at t-1, X(t-1) 115 Pt_1 = M(1:end-1); 116 117 % Discretization for daily prices 118 dt = 1/365; 119 120 % PDF for discretized model 121 mrjpdf = @(Pt, phi, sigmaSq) ... 122 exp((-(Pt-phi.*Pt_1).^2)/(2.*sigmaSq)).* ... 123 (1/sqrt(2.*pi.*sigmaSq)); 124 125 % Constraints: 126 lb = [-Inf 0]; 127 ub = [1 Inf]; 128 129 % Initial values 130 x0 = [0 var(M)]; 131 132 % Maximum likelihood estimation 133 params_Msa = mle(Pt,’pdf’,mrjpdf,’start’,x0,’lowerbound’,lb,’upperbound’,ub,... 134 ’optimfun’,’fmincon’); 135 136 % Calibrated parameters 137 138 kappa_Msa = params_Msa(1) 139 140 sigma_Msa = sqrt(params_Msa(2)) 141 142 143 --------------------------raintemp-------------------------- 144 145 %clear; clc; close all; 146 147 Z = xlsread(’kiambu_rain_temp_data.xlsx’, ’Sheet1’); 148 149 DateTime = Z(:,1); 150 Rain = Z(:,2); 151 Temp = Z(:,3); 152 153 MATLABDates = x2mdate(DateTime,0,’datetime’) 154 155 save(’MATLABDates.mat’) 156 save(’Rain.mat’) 157 save(’Temp.mat’) 158 159 160 --------------------------kirinyaga_temp-------------------------- 161 162 %clear; clc; close all; 163 164 load(’Temp.mat’) 165 x = length(Temp); 166 t = (1:1:x)’; 167 168 init = [20, 1, 1, 10, 1, 10]; 169 options = optimoptions(@lsqnonlin,’Algorithm’,’trust-region-reflective’); 170 modelfun = @(b) b(1) + b(2)*t + b(3)*cos(2*pi*(t - b(4))/365)... 32 171 + b(5)*(1 + sin(2*pi*(t - b(6))/365)).*sin(2*pi*t/365) - Temp; 172 x = lsqnonlin(modelfun,init,[],[],options) 173 Alpha = TruncatedF(x(1), x(2), x(3), x(4), x(5), x(6), t); 174 175 E=Temp-Alpha 176 177 E 178 179 [B,BINT,V2] = regress(E(2:end),E(1:end-1)) 180 181 182 --------------------------TruncatedFTemp-------------------------- 183 184 function Alpha = TruncatedF(a0,a1,a2,a3,a4,a5,t) 185 % Seasonal function with trend, the truncated Fourier transform 186 %Gamma = TruncatedF(x(1), x(2), x(3), x(4), x(5), x(6), t); 187 Alpha = a0 + a1*t + a2*cos(2*pi*(t-a3)/365)+ a4*(1+sin(2*pi*(t - a5)/365)).*sin(2*pi*t/365); 188 end 189 190 191 --------------------------antitheticMC-------------------------- 192 193 function [call, put ] = antitheticMCtrial (K, TiN, TiM, S, kappa_Nai, kappa_Msa, 194 sigma_Nai, sigma_Msa, h, T, rho) 195 K = 500; 196 N = 10000000; 197 r=0.07 198 h=1/N; 199 T=1; 200 rho=0.6060 201 202 TiN=6000; 203 TiM=5000; 204 205 a0 = 4.9743; 206 a1 = -0.0001; 207 a2 = -0.8866; 208 a3 = -2.4918; 209 a4 = 0.7568; 210 a5 = -0.0273; 211 %n=50; 212 213 for i=1:N 214 S(i,:)=a0 + a1*i + a2*cos(2*pi*(i-a3)/365)+ a4*(1+sin(2*pi*(i - a5)/365)).*sin(2*pi*i/365); 215 %a0+a1*(i,:)+a2*cos*(2*pi S(i,:).*(1+r*h+sigma*dW(i,:)); 216 end 217 218 % We generate two standard Gaussian vectors 219 U = randn(N,1); V = randn(N,1); 220 % We generate the antithetic random vectors 221 nU = -U ; nV = -V; 222 223 224 225 % We generate final prices for both assets from both random vectors 226 pTi_1 = TiN+ S(i,:)-S(i)-kappa_Nai*(TiN-S(i))*T+sigma_Nai.*sqrt(T).*U 227 nTi_1 = TiN+ S(i,:)-S(i)-kappa_Nai*(TiN-S(i))*T+sigma_Nai.*sqrt(T).*nU 228 pTi_2 = TiM+ S(i,:)-S(i)-kappa_Msa*(TiM-S(i))*T+sigma_Msa.* 33 229 (sqrt(T).*V.*rho+sqrt(1-rho).*sqrt(T).*V) 230 nTi_2 = TiM+ S(i,:)-S(i)-kappa_Msa*(TiM-S(i))*T+sigma_Msa.* 231 (sqrt(T).*nV.*rho+sqrt(1-rho).*sqrt(T).*nV) 232 233 % We compute the payoff vector for the call for both random vectors 234 resCp = max((pTi_1 - pTi_2) - K, 0); 235 resCn = max((nTi_1 - nTi_2) - K, 0); 236 % We compute the payoff vector for the put for both random vectors 237 resPp = max(K - (pTi_1 - pTi_2), 0); 238 resPn = max(K - (nTi_1 - nTi_2), 0); 239 240 % We compute the average between normal and antithetic payoffs 241 resC = 0.5*(resCp + resCn); 242 resP = 0.5*(resPp + resPn); 243 244 % We finally discount the average of the payoff 245 call = exp(-r*T) * mean(resC); 246 put = exp(-r*T) * mean(resP); 247 end 248 249 250 --------------------------stderrorMC-------------------------- 251 252 sqrt((sigma_Nai)^2/10) 253 sqrt((sigma_Nai)^2/100) 254 sqrt((sigma_Nai)^2/1000) 255 sqrt((sigma_Nai)^2/10000) 256 sqrt((sigma_Nai)^2/100000) 257 sqrt((sigma_Nai)^2/1000000) 258 sqrt((sigma_Nai)^2/10000000) 259 260 sqrt((sigma_Msa^2)/10) 261 sqrt((sigma_Msa^2)/100) 262 sqrt((sigma_Msa^2)/1000) 263 sqrt((sigma_Msa^2)/10000) 264 sqrt((sigma_Msa^2)/100000) 265 sqrt((sigma_Msa^2)/1000000) 266 sqrt((sigma_Msa^2)/10000000) 34