-
Notifications
You must be signed in to change notification settings - Fork 0
/
gas_dejong.py
190 lines (121 loc) · 8 KB
/
gas_dejong.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
# -*- coding: utf-8 -*-
"""
Created on Fri Mar 9 10:56:16 2018
@author: fabdellah
"""
import numpy as np
from numpy import random
def payoff1(S,dv,injection_cost,withdrawal_cost):
if dv > 0:
h= -injection_cost*dv
elif dv == 0:
h = 0
else:
h = -withdrawal_cost*dv
return h
def payoff(S,dv,injection_cost,withdrawal_cost):
return -withdrawal_cost*dv
def penalty(S,v):
#specify the penalty function otherwise return 0
return 0
class gas_storage(object):
"""
S0 : initial spot price
T : time to maturity
steps : number of discrete times (delta_t = T/steps)
r : riskless discount rate (constant)
sigma : volatility of returns
"""
def __init__(self, S0, V0, V_end, T, steps,M, r, sigma, nbr_simulations,vMax,vMin , max_rate, min_rate , injection_cost, withdrawal_cost ):
self.S0 = S0 # Parameters for the spot price
self.T = T
self.r = r
self.sigma = sigma
self.nbr_simulations = nbr_simulations
self.V0 = V0
self.V_end = V_end
self.vMax = vMax # Parameters for facility
self.vMin = vMin
self.max_rate = max_rate
self.min_rate = min_rate
self.injection_cost =injection_cost # Parameters for the payoff function
self.withdrawal_cost = withdrawal_cost
if S0 < 0 or T <= 0 or r < 0 or sigma < 0 or nbr_simulations < 0 :
raise ValueError('Error: Negative inputs not allowed')
self.steps = steps
self.M = M
self.alpha = self.vMax/(self.M-1)
self.delta_t = self.T / float(self.steps)
self.discount = np.exp(-self.r * self.delta_t)
def simulated_price_matrix(self, seed = 1):
""" Returns Monte Carlo simulated prices (matrix)
rows: time
columns: price-path simulation """
np.random.seed(seed)
simulated_price_matrix = np.zeros((self.steps + 2, self.nbr_simulations), dtype=np.float64)
simulated_price_matrix[0,:] = self.S0
for t in range(1, self.steps + 2):
brownian = np.random.standard_normal( int(self.nbr_simulations / 2))
brownian = np.concatenate((brownian, -brownian))
simulated_price_matrix[t, :] = (simulated_price_matrix[t - 1, :]
* np.exp((self.r - self.sigma ** 2 / 2.) * self.delta_t
+ self.sigma * brownian * np.sqrt(self.delta_t)))
#needs to be specified according to the corresponding 2-factor model
return simulated_price_matrix
def contract_value(self):
value_matrix = np.zeros((self.simulated_price_matrix().shape[0],self.simulated_price_matrix().shape[1],self.M)) # time, path , volume level
acc_cashflows = np.zeros_like(value_matrix)
decision_rule = np.zeros_like(value_matrix)
volume_level = np.zeros_like(self.simulated_price_matrix())
decision_rule_avg = np.zeros((self.simulated_price_matrix().shape[0],self.M))
volume_level_avg = np.zeros(self.simulated_price_matrix().shape[0])
acc_cashflows_avg = np.zeros((self.simulated_price_matrix().shape[0],self.M))
volume_level_avg[0] = self.V0
#volume_level_avg[-1] = self.V_end
volume_level[0, : ] = self.V0*np.ones(self.nbr_simulations)
#volume_level[-1,: ] = self.V_end*np.ones(self.nbr_simulations)
value_matrix[-1,: ,:] = penalty(self.simulated_price_matrix()[-1, :],volume_level[-1,:])
acc_cashflows[-1,:,:] = penalty(self.simulated_price_matrix()[-1, :],volume_level[-1,:])
from scipy.optimize import minimize
for t in range(self.steps , 0 , -1):
print ('-----------')
print ('Time: %5.3f, Spot price: %5.1f ' % (t, self.simulated_price_matrix()[t, 1]))
for m in range(1,self.M):
volume_level[t+1,:] = (m-1)*self.alpha
regression = np.polyfit(self.simulated_price_matrix()[t, :], acc_cashflows[t+1, :, m-1] * self.discount, 2)
continuation_value = np.polyval(regression, self.simulated_price_matrix()[t, :])
for b in range(self.nbr_simulations):
f = lambda x: -1*( payoff(self.simulated_price_matrix()[t, b],
x ,self.injection_cost,self.withdrawal_cost ) + continuation_value[b] )
cons = ({'type': 'ineq', 'fun': lambda x: (volume_level[t+1,b] - x - self.vMin) },
{'type': 'ineq', 'fun': lambda x: (self.vMax - volume_level[t+1,b] + x) },
{'type': 'ineq', 'fun': lambda x: (self.max_rate - x) },
{'type': 'ineq', 'fun': lambda x: (x - self.min_rate) },
{'type': 'ineq', 'fun': lambda x: (self.max_rate*t - volume_level[t+1,b] + x ) })
res = minimize(f, random.rand(1), constraints=cons)
decision_rule[t,b,m-1] = res.x
#if (decision_rule[t,b,m] != 0):
#volume_level[t,b] = volume_level[t+1,b]- (np.sign(decision_rule[t,1,m]))*((abs(decision_rule[t,1,m])//self.alpha)+1)*self.alpha
#else:
#volume_level[t,b] = volume_level[t+1,b]
#print ('Time: %5.2f, Spot price: %5.1f , Decision rule (inject/withdraw): %5.1f , Volume level: %5.3f, Acc_cashflows: %5.2f '
# % (t, self.simulated_price_matrix()[t, b] , decision_rule[t,b,m], volume_level[t,b], acc_cashflows[t,b,m] ))
acc_cashflows[t,:,m-1] = payoff(self.simulated_price_matrix()[t, :],
decision_rule[t,:,m-1] , self.injection_cost, self.withdrawal_cost) + acc_cashflows[t+1,:,m-1]*self.discount
decision_rule_avg[t,m-1] = np.sum(decision_rule[t,:,m-1])/self.nbr_simulations
volume_level_avg[t] = np.sum(volume_level[t,:])/self.nbr_simulations
acc_cashflows_avg[t,m-1] = np.sum(acc_cashflows[t,:,m-1])/self.nbr_simulations
print (' Decision rule (inject/withdraw): %5.1f , Acc_cashflows: %5.2f , payoff: %5.2f '
% (decision_rule[t,1,m], acc_cashflows[t,1,m-1], payoff(self.simulated_price_matrix()[t, 1],decision_rule[t,1,m-1] , self.injection_cost, self.withdrawal_cost)))
#print('sign: %5.3f , deux: %5.3f , v-add: %5.3f ' % ( np.sign(decision_rule[t,1,m]) , (abs(decision_rule[t,1,m])//self.alpha)+1 , (np.sign(decision_rule[t,1,m]))*((abs(decision_rule[t,1,m])//self.alpha)+1)*self.alpha ))
return acc_cashflows[1,:,:] * self.discount # at time 0
def price(self):
return round((np.sum(self.contract_value(),axis=0) / float(self.nbr_simulations))[0] , 2)
# gas_storage(S0, V0,V_end, T, steps, M, r, sigma, nbr_simulations, vMax, vMin , max_rate , min_rate , injection_cost, withdrawal_cost )
facility1 = gas_storage(25, 0 , 0, 1, 12, 101, 0.06, 0.4, 10, 5000000, 0 , 34000, -53000 ,0.1, 0.1)
from time import time
t0 = time()
print ('Price: ' ,facility1.price() )
t1 = time();
d1 = t1 - t0
print ("Duration in Seconds %6.3f" % d1)