Coverage for rivapy/models/residual_demand_fwd_model.py: 45%
336 statements
« prev ^ index » next coverage.py v7.8.2, created at 2025-06-05 14:27 +0000
1import numpy as np
2import pandas as pd
3import datetime as dt
4import matplotlib.pyplot as plt
5import abc
6from typing import Union, Callable, List, Tuple, Dict, Protocol, Set
7import scipy
8from scipy.special import comb
9from rivapy.tools.interfaces import FactoryObject
10from rivapy.models.factory import create as _create
11from rivapy.models.ornstein_uhlenbeck import OrnsteinUhlenbeck
12from rivapy.models.base_model import BaseFwdModel
14def _logit(x):
15 return np.log(x/(1-x))
17def _inv_logit(x):
18 return 1.0/(1+np.exp(-x))
20class ForwardSimulationResult(abc.ABC):
21 @abc.abstractmethod
22 def n_forwards(self)->float:
23 pass
25 @abc.abstractmethod
26 def udls(self)->Set[str]:
27 pass
29 def keys(self)->List[str]:
30 result = set()
31 for udl in self.udls():
32 for i in range(self.n_forwards()):
33 result.add(BaseFwdModel.get_key(udl, i))
34 return result
36 @abc.abstractmethod
37 def get(self, key: str)->np.ndarray:
38 pass
40class WindPowerForecastModelParameter(FactoryObject):
41 def __init__(self, n_call_strikes:int = 20, min_strike: float=-5.0, max_strike: float=5.0):
42 """Parameter class for the wind power forecast model.
44 Args:
45 n_call_strikes (int, optional): Number of calls used to approximate the expectation of inverse logit applied to Ornsetin-Uhlenbeck. Defaults to 20.
46 min_strike (float, optional): Minimal strike used to compute the approximation of inverse logit applied to Ornsetin-Uhlenbeck. Defaults to -5.0.
47 max_strike (float, optional): Maximal strike used to compute the approximation of inverse logit applied to Ornsetin-Uhlenbeck. Defaults to 5.0.
48 """
49 self.n_call_strikes = n_call_strikes
50 self.min_strike = min_strike
51 self.max_strike = max_strike
53 def _to_dict(self) -> dict:
54 return {'n_call_strikes': self.n_call_strikes,
55 'min_strike': self.min_strike,
56 'max_strike': self.max_strike}
58class WindPowerForecastModel(BaseFwdModel):
61 class ForwardSimulationResult(ForwardSimulationResult):
62 def __init__(self, paths:np.ndarray,
63 wind_forecast_model,
64 timegrid,
65 expiries,
66 initial_forecasts):
67 self._paths = paths
68 self._timegrid = timegrid
69 self._model = wind_forecast_model
70 self._result = None
71 self.expiries = expiries
72 self.initial_forecasts = initial_forecasts
73 self._ou_additive_forward_corrections = wind_forecast_model._compute_ou_additive_forward_correction(self.expiries, self.initial_forecasts)
75 def n_forwards(self)->float:
76 return len(self.expiries)
78 def udls(self)->Set[str]:
79 return self._model.udls()
81 def get(self, key: str,
82 forecast_timepoints: List[int]=None)->np.ndarray:
83 expiry = BaseFwdModel.get_expiry_from_key(key)
84 return self._get(expiry, forecast_timepoints)
86 def _get(self, expiry: int, forecast_timepoints: List[int])->np.ndarray:
87 result = np.empty(self._paths.shape)
88 for i in range(self._timegrid.shape[0]):
89 #result[i,:] = self._model.get_forward(self._paths[i,:], self.expiries[expiry]-self._timegrid[i], self._ou_additive_forward_corrections[expiry])
90 self._model._compute_expectation_inv_logit(result[i,:],
91 self._paths[i,:],
92 self.expiries[expiry]-self._timegrid[i],
93 self._ou_additive_forward_corrections[expiry])
94 if forecast_timepoints is not None:
95 ftp_prev = 1
96 for ftp in forecast_timepoints:
97 for i in range(ftp_prev, ftp):
98 if abs(self._timegrid[i]-self.expiries[expiry])>1e-6:
99 result[i,:] = result[i-1,:]
100 ftp_prev = ftp+1
101 for i in range(ftp_prev, self._timegrid.shape[0]):
102 if abs(self._timegrid[i]-self.expiries[expiry])>1e-6:
103 result[i,:] = result[i-1,:]
104 return result
106 def __init__(self, region: str,
107 speed_of_mean_reversion: float,
108 volatility: float,
109 params: WindPowerForecastModelParameter=None
110 ):
111 """Simple model to simulate forecasts of wind efficiencies (power production by wind as percentage of total wind capacity) based on the Ornstein-Uhlenbeck process.
113 The base model used within this model is the Ornstein-Uhlenbeck process and uses internally the class :class:`rivapy.models.OrnsteinUhlenbeck` with a mean level of zero to simulate the forecast.
114 Here, the forecast of wind efficiency :math:`w_{t,T}` at time :math:`T` is computed as the expected value of the standard logistic function (sigmoid) applied to the Ornstein-Uhlenbeck process conditioned on current simulated value,
115 i.e.
117 .. math:: w_{t, T} = \\frac{1}{1+e^{-(X_{t,T}+\\nu(T))}}
119 where :math:`X_{t,T}:=E[X_T\mid X_t]` with
121 .. math:: dX = -\\lambda X dt + \\sigma dW_t, X(0) = 0
123 and :math:`\\nu(T)` is a correction term to ensure that the initial value of :math:`w_{t,T}` is equal to
124 a given initial forecast :math:`\\bar{w}_{0,T}`. It is computed by
126 .. math:: \\nu(T) := \log\\frac{\\bar{w}_{0,T}}{(1.0-\\bar{w}_{0,T})}.
128 The sigmoid function is applied to ensure that the forecast is between 0 and 1 (as percentage of total wind capacity)
130 Args:
131 speed_of_mean_reversion (float): The speed of mean reversion of the underlying Ornstein-Uhlenbeck process (:class:`rivapy.models.OrnsteinUhlenbeck`).
132 volatility (float): The volatility of the underlying Ornstein-Uhlenbeck process (:class:`rivapy.models.OrnsteinUhlenbeck`).
133 region (str): The name of the respective wind region.
134 """
135 #if len(expiries) != len(initial_forecasts):
136 # raise Exception('Number of forward expiries does not equal number of initial forecasts. Each forward expiry needs an own forecast.')
137 self.ou = OrnsteinUhlenbeck(speed_of_mean_reversion, volatility, mean_reversion_level=0.0)
138 self.region = region
139 if params is None:
140 self.params = WindPowerForecastModelParameter()
141 else:
142 self.params = params
143 self.call_strikes, self.call_weights = WindPowerForecastModel._compute_strikes_weights(self.params.n_call_strikes,
144 self.params.min_strike,
145 self.params.max_strike)
146 def _compute_expectation_inv_logit(self, result: np.ndarray, spots: np.ndarray,
147 ttm: float, strike_offset: float=0.0):
148 """Compute the expectation of the inverse logit applied to the ou model.
150 Args:
151 result (np.ndarray): Array to store the result
152 spots (np.ndarray): Array of spots
153 ttm (float): Time to maturity
154 strike_offset (float, optional): Offset to the strike (needed to fit the forecast at initial time). Defaults to 0.0.
155 """
156 strikes = self.call_strikes
157 ref_spots = np.linspace(spots.min(), spots.max(), num=self.params.n_call_strikes, endpoint=True)
158 prices = np.empty((strikes.shape[0], self.params.n_call_strikes,))
159 for i in range(strikes.shape[0]):
160 prices[i,:] = self.ou.compute_call_price(ref_spots, strikes[i]-strike_offset, ttm)
162 result[:] = self.call_weights[0]*np.interp(spots, ref_spots, prices[0,:])
163 for i in range(1, prices.shape[0]):
164 result += self.call_weights[i]*np.interp(spots, ref_spots, prices[i,:])
166 def _compute_ou_additive_forward_correction(self, expiries, initial_forecasts):
167 # we determine the additive correction term such that the expected value
168 # of the inverse logit is equal to the initial forecast
169 def error(i):
170 def _error(correction):
171 tmp = 0
172 for j in range(self.call_strikes.shape[0]):
173 tmp += self.call_weights[j]*self.ou.compute_call_price(0.0, self.call_strikes[j]-correction, expiries[i])
174 return tmp-initial_forecasts[i]
175 return _error
177 result = np.empty((len(expiries)))
178 for i in range(len(expiries)):
179 result[i] = scipy.optimize.brentq(error(i), -15.0, 15.0, xtol=1e-6)
180 return result
182 def _to_dict(self)->dict:
183 return {'speed_of_mean_reversion': self.ou.speed_of_mean_reversion,
184 'volatility': self.ou.volatility,
185 'region': self.region}
187 @staticmethod
188 def _compute_strikes_weights(n_strikes: int=20, min_strike: float=-5.0, max_strike: float=5.0):
189 strikes = np.linspace(start=min_strike, stop=max_strike, num=n_strikes)
190 f = _inv_logit(strikes)
191 weights = np.zeros((strikes.shape[0]+1))
192 h = strikes[1]-strikes[0]
193 butterfly_weights = [1, -2, 1]/h
194 for s in range(1,strikes.shape[0]):
195 weights[s-1] += f[s]*butterfly_weights[0]
196 weights[s] += f[s]*butterfly_weights[1]
197 weights[s+1] += f[s]*butterfly_weights[2]
198 weights[-2] -= 0.5* butterfly_weights[1]*f[s]
199 return strikes, weights[:-1]
202 @staticmethod
203 def eval_call_functions(strikes, weights, x):
204 approx = np.copy(weights[0]*np.maximum(x-strikes[0],0.0))
205 for i in range(1,weights.shape[0]):
206 approx += weights[i]*np.maximum(x-strikes[i],0.0)
207 return approx
209 def get_forward(self, paths, ttm, ou_additive_forward_correction: float):
210 expected_ou = self.ou.compute_expected_value(paths, ttm)#+correction
211 result = _inv_logit(expected_ou + ou_additive_forward_correction)
212 #expected_ou = self.ou.compute_expected_value(paths[1], ttm)#+correction
213 #result += (1.0-ttm)*_inv_logit(expected_ou + ou_additive_forward_correction)
214 return result
216 def rnd_shape(self, n_sims: int, n_timesteps: int)->tuple:
217 return (n_timesteps-1, n_sims)
219 def simulate(self, timegrid,
220 rnd: np.ndarray,
221 expiries: List[float],
222 initial_forecasts: List[float],
223 startvalue=0.0)->ForwardSimulationResult:
224 #paths = np.empty((2,timegrid.shape[0], rnd.shape[2]))
225 paths = self.ou.simulate(timegrid, startvalue, rnd)
226 #paths[1,:,:] = self.ou.simulate(timegrid, startvalue, rnd[1])
227 return WindPowerForecastModel.ForwardSimulationResult(paths, self, timegrid, expiries, initial_forecasts)
229 def udls(self)->Set[str]:
230 return set([self.region])
232class MultiRegionWindForecastModel(BaseFwdModel):
233 class ForwardSimulationResult(ForwardSimulationResult):
234 def __init__(self, model, regions_result):
235 self._results = regions_result
236 self._model = model
238 def n_forwards(self)->float:
239 for v in self._results.values():
240 return v.n_forwards()
242 def udls(self)->Set[str]:
243 return self._model.udls()
245 def get(self, key: str, forecast_timepoints: List[int]=None)->np.ndarray:
246 udl = BaseFwdModel.get_udl_from_key(key)
247 if udl == self._model.name:
248 expiry = BaseFwdModel.get_expiry_from_key(key)
249 return self._get(expiry, forecast_timepoints)
250 for v in self._results.values():
251 if udl in v.udls():
252 return v.get(key, forecast_timepoints)
253 raise Exception('Cannot find key '+ key)
255 def _get(self, expiry: int, forecast_timepoints: List[int])->np.ndarray:
256 result = None
257 for udl in self.udls():
258 if udl == self._model.name:
259 continue
260 if result is None:
261 result = self._model.region_relative_capacity(udl)*np.copy(self.get(BaseFwdModel.get_key(udl, expiry), forecast_timepoints))
262 else:
263 result += self._model.region_relative_capacity(udl)*np.copy(self.get(BaseFwdModel.get_key(udl, expiry), forecast_timepoints))
264 return result
266 class Region(FactoryObject):
267 def __init__(self, model: WindPowerForecastModel, capacity: float, rnd_weights: List[float]):
268 self.model = _create(model)
269 self.capacity = capacity
270 self.rnd_weights = rnd_weights
272 def name(self):
273 return self.model.region
275 def n_random(self):
276 return len(self.rnd_weights)
278 def udls(self)->Set[str]:
279 return self.model.udls()
281 def _to_dict(self)->dict:
282 return {'model': self.model.to_dict(), 'capacity': self.capacity, 'rnd_weights': self.rnd_weights}
284 def __init__(self, name: str, region_forecast_models: List[Region]):
285 """Simple model to simulate power forecasts for more then one wind region and a total efficiency over all regions.
287 This model relates the wind models for the different regions within the simulation by just using a linear sum of normal random variates, i.e.
288 for a region :math:`r` we have weights :math:`w_{r,i}`, :math:`1\leq i\leq N`, so that within the simulation based on :math:`n` different normal random
289 variates :math:`X_i` we compute the random variate for the region by
291 .. math:: X_{r,i} = \\frac{\sum_i w_{r,i} X_i}{\sqrt{\sum_i w_{r,i}^2}}.
293 Args:
294 name (str): Name of the overall region.
295 region_forecast_models (List[Region]): List of regions that will be simulated.
297 Examples:
298 >>> wind_onshore = WindPowerForecastModel(region='Onshore', speed_of_mean_reversion=0.1, volatility=4.80)
299 >>> wind_offshore = WindPowerForecastModel(region='Offshore', speed_of_mean_reversion=0.5, volatility=4.80)
300 >>> regions = [ MultiRegionWindForecastModel.Region(
301 wind_onshore,
302 capacity=1000.0,
303 rnd_weights=[0.8,0.2]
304 ),
305 MultiRegionWindForecastModel.Region(
306 wind_offshore,
307 capacity=100.0,
308 rnd_weights=[0.2,0.8]
309 )
310 ]
311 >>> wind = MultiRegionWindForecastModel('Wind_Germany', regions)
312 >>> # after model setup we simulate forecasts
313 >>> days = 10
314 >>> timegrid = np.linspace(0.0, days*1.0/365.0, days*24)
315 >>> forward_expiries = [timegrid[-1] + i/(365.0*24.0) for i in range(4)]
316 >>> rnd = np.random.normal(size=wind.rnd_shape(n_sims, timegrid.shape[0]))
317 >>> results = wind.simulate(timegrid, rnd, expiries=forward_expiries,
318 initial_forecasts={'Onshore': [0.8, 0.7,0.6,0.5],
319 'Offshore': [0.6, 0.6, 0.5, 0.5]}
320 )
321 """
322 self.name = name
323 if len(region_forecast_models)==0:
324 raise Exception('Empty list of models is not allowed')
325 self._region_forecast_models = [ _create(r) for r in region_forecast_models]
326 n_rnd_ref_model = self._region_forecast_models[0].n_random()
327 for i in range(1, len(self._region_forecast_models)):
328 if self._region_forecast_models[i].n_random() != n_rnd_ref_model:
329 raise Exception('All regions must have the same number of random variables.')
333 def _to_dict(self)->dict:
334 return {'name': self.name,
335 'region_forecast_models':[v.to_dict() for v in self._region_forecast_models]
336 }
338 def rnd_shape(self, n_sims: int, n_timesteps: int)->tuple:
339 return (self._region_forecast_models[0].rnd_weights, n_timesteps-1, n_sims)
341 def total_capacity(self):
342 result = 0.0
343 for r in self._region_forecast_models:
344 result += r.capacity
345 return result
347 def region_relative_capacity(self, region: str):
348 for r in self._region_forecast_models:
349 if r.name() == region:
350 return r.capacity/self.total_capacity()
351 raise Exception('Model does not contain a region with name ' + region)
353 def region_names(self)->List[str]:
354 return [r.name() for r in self._region_forecast_models]
356 def simulate(self, timegrid, rnd: np.ndarray, expiries: List[float],
358 initial_forecasts: Dict[str,List[float]], startvalue=0.0):
359 results = {}
360 for region in self._region_forecast_models:
361 rnd_ = region.rnd_weights[0]*rnd[0,:,:]
362 for i in range(1, region.n_random()):
363 rnd_ += region.rnd_weights[i]*rnd[i,:,:]
364 results[region.name()] = region.model.simulate(timegrid, rnd_, expiries, initial_forecasts[region.name()], startvalue)
365 return MultiRegionWindForecastModel.ForwardSimulationResult(self, results)
367 def rnd_shape(self, n_sims: int, n_timesteps: int)->tuple:
368 return (len(self._region_forecast_models), n_timesteps-1, n_sims)
370 def udls(self)->Set[str]:
371 result = set([self.name])
372 for v in self._region_forecast_models:
373 result.update(v.udls())
374 return result
376class LinearDemandForwardModel(BaseFwdModel):
377 class ForwardSimulationResult(ForwardSimulationResult):
378 def __init__(self, model, highest_price, wind_results, additive_correction):
379 self._model = model
380 self._highest_price = highest_price
381 self._wind = wind_results
382 self._additive_correction = additive_correction
384 def n_forwards(self)->float:
385 return self._wind.n_forwards()
387 def udls(self)->Set[str]:
388 return self._model.udls()
390 def get(self, key: str, forecast_timepoints: List[int]=None)->np.ndarray:
391 udl = BaseFwdModel.get_udl_from_key(key)
392 if udl == self._model.power_name:
393 expiry = BaseFwdModel.get_expiry_from_key(key)
394 return self._get(expiry, forecast_timepoints)
395 else:
396 return self._wind.get(key, forecast_timepoints)
398 def _get(self, expiry: int, forecast_timepoints: List[int])->np.ndarray:
399 total_produced = self._wind.get(BaseFwdModel.get_key(self._wind._model.name, expiry), forecast_timepoints)
400 power_fwd = np.empty((total_produced.shape[0], total_produced.shape[1]))
401 for i in range(total_produced.shape[0]):
402 power_fwd[i,:] = (1.0-total_produced[i,:] )*(self._highest_price[i,:]+self._additive_correction[expiry])
403 return power_fwd
405 def __init__(self, wind_power_forecast: MultiRegionWindForecastModel,
406 x_volatility: float,
407 x_mean_reversion_speed: float,
408 power_name:str = None):
409 self.wind_power_forecast = _create(wind_power_forecast)
410 self.highest_price_ou_model: OrnsteinUhlenbeck = OrnsteinUhlenbeck(x_mean_reversion_speed, x_volatility, 0.0)
411 if power_name is not None:
412 self.power_name = power_name
413 else:
414 self.power_name = 'POWER'
415 #self.region_to_capacity = region_to_capacity
417 def _to_dict(self)->dict:
418 return {'wind_power_forecast': self.wind_power_forecast.to_dict(),
419 'x_volatility': self.highest_price_ou_model.volatility,
420 'x_mean_reversion_speed': self.highest_price_ou_model.speed_of_mean_reversion,
421 'power_name': self.power_name,
422 }
424 def rnd_shape(self, n_sims: int, n_timesteps: int)->tuple:
425 rnd_shape = self.wind_power_forecast.rnd_shape(n_sims, n_timesteps)
426 if len(rnd_shape) == 3:
427 return (rnd_shape[0]+1, rnd_shape[1], rnd_shape[2])
428 return (2, rnd_shape[0], rnd_shape[1])
430 def get_technology(self)->str:
431 """Return name of the technology modeled.
433 Returns:
434 str: Name of instrument.
435 """
436 return self.wind_power_forecast.region
438 def _compute_additive_correction(self, expiries: List[float],
439 power_fwd_prices: List[float],
440 initial_forecasts: Dict[str,List[float]], startvalue=0.0 ):
441 result = np.empty((len(expiries)))
442 for i in range(len(expiries)):
443 tmp = 0
444 for k,v in initial_forecasts.items():
445 tmp += v[i]*self.wind_power_forecast.region_relative_capacity(k)
446 if abs(1.0-tmp) < 1e-5:
447 raise Exception('Initial forecasts sum to 1.0. Cannot compute additive correction.')
448 result[i] = power_fwd_prices[i]/(1.0-tmp) - self.highest_price_ou_model.compute_expected_value(startvalue, expiries[i])
449 return result
451 def simulate(self, timegrid: np.ndarray,
452 rnd: np.ndarray,
453 expiries: List[float],
454 power_fwd_prices: List[float],
455 initial_forecasts: Dict[str, List[float]]):
456 additive_correction = self._compute_additive_correction(expiries, power_fwd_prices, initial_forecasts, startvalue=0.0)
457 highest_prices = self.highest_price_ou_model.simulate(timegrid, 0.0, rnd[0,:])
458 simulated_wind = self.wind_power_forecast.simulate(timegrid, rnd[1:,:], expiries, initial_forecasts)
459 return LinearDemandForwardModel.ForwardSimulationResult(self, highest_prices, simulated_wind, additive_correction)
461 def udls(self)->Set[str]:
462 result = self.wind_power_forecast.udls()
463 result.add(self.power_name)
464 return result
469class ResidualDemandForwardModel(BaseFwdModel):
470 class ForwardSimulationResult(ForwardSimulationResult):
471 def __init__(self, model, highest_price, wind_results):
472 self._model = model
473 self._highest_price = highest_price
474 self._wind = wind_results
476 def n_forwards(self)->float:
477 return self._wind.n_forwards()
479 def udls(self)->Set[str]:
480 return self._model.udls()
482 def get(self, key: str, forecast_timepoints: List[int])->np.ndarray:
483 udl = BaseFwdModel.get_udl_from_key(key)
484 if udl == self._model.power_name:
485 expiry = BaseFwdModel.get_expiry_from_key(key)
486 return self._get(expiry, forecast_timepoints)
487 else:
488 return self._wind.get(key, forecast_timepoints)
490 def _get(self, expiry: int, forecast_timepoints: List[int])->np.ndarray:
491 total_produced = self._wind.get(BaseFwdModel.get_key(self._wind._model.name, expiry), forecast_timepoints)
492 power_fwd = np.empty((total_produced.shape[0], total_produced.shape[1]))
493 for i in range(total_produced.shape[0]):
494 power_fwd[i,:] = self._model.supply_curve(1.0-total_produced[i,:] )*self._highest_price[i,:]
495 return power_fwd
497 def __init__(self, wind_power_forecast,
498 highest_price_ou_model,
499 supply_curve: Callable[[float], float],
500 max_price: float,
501 power_name:str = None):
502 #print(wind_power_forecast)
503 self.wind_power_forecast = _create(wind_power_forecast)
504 self.highest_price_ou_model = _create(highest_price_ou_model)
505 self.supply_curve = _create(supply_curve)
506 self.max_price = max_price
507 if power_name is not None:
508 self.power_name = power_name
509 else:
510 self.power_name = 'POWER'
511 #self.region_to_capacity = region_to_capacity
513 def _to_dict(self)->dict:
514 return {'wind_power_forecast': self.wind_power_forecast.to_dict(),
515 'supply_curve': self.supply_curve.to_dict(),
516 'highest_price_ou_model': self.highest_price_ou_model.to_dict(),
517 'max_price': self.max_price,
518 'power_name': self.power_name,
519 #'region_to_capacity': self.region_to_capacity
520 }
522 def rnd_shape(self, n_sims: int, n_timesteps: int)->tuple:
523 rnd_shape = self.wind_power_forecast.rnd_shape(n_sims, n_timesteps)
524 if len(rnd_shape) == 3:
525 return (rnd_shape[0]+1, rnd_shape[1], rnd_shape[2])
526 return (2, rnd_shape[0], rnd_shape[1])
528 def get_technology(self)->str:
529 """Return name of the technology modeled.
531 Returns:
532 str: Name of instrument.
533 """
534 return self.wind_power_forecast.region
536 def simulate(self, timegrid: np.ndarray,
537 rnd: np.ndarray,
538 expiries: List[float],
539 initial_forecasts: Dict[str, List[float]]):
540 highest_prices = self.highest_price_ou_model.simulate(timegrid, 1.0, rnd[0,:])*self.max_price
541 simulated_wind = self.wind_power_forecast.simulate(timegrid, rnd[1:,:], expiries, initial_forecasts)
542 return ResidualDemandForwardModel.ForwardSimulationResult(self, highest_prices, simulated_wind)
544 def udls(self)->Set[str]:
545 result = self.wind_power_forecast.udls()
546 result.add(self.power_name)
547 return result
551if __name__=='__main__':
552 params = WindPowerForecastModelParameter(n_call_strikes=40, min_strike=-7.0, max_strike=7.0)
553 # setup
554 wind_region_model = {}
555 vols = [1.0]
556 mean_reversion_speed = [0.5]
557 capacities = [10_000.0]
558 rnd_weights = [ [1.0]
559 ]
560 np.random.seed(42)
561 regions = []
562 for i in range(len(vols)):
563 model = WindPowerForecastModel(region='Region_' + str(i),
564 speed_of_mean_reversion=mean_reversion_speed[i],
565 volatility=vols[i], params=params)
566 regions.append(MultiRegionWindForecastModel.Region(
567 model,
568 capacity=capacities[i],
569 rnd_weights=rnd_weights[i]
570 ) )
571 wind = MultiRegionWindForecastModel('Wind_Germany', regions)
572 highest_price = OrnsteinUhlenbeck(speed_of_mean_reversion=1.0, volatility=0.5, mean_reversion_level=0.0)
573 model = LinearDemandForwardModel(wind_power_forecast=wind,
574 highest_price_ou_model= highest_price,
575 power_name= 'Power_Germany')
576 timegrid = np.linspace(0.0,1.0, 365)
577 np.random.seed(42)
578 rnd = np.random.normal(size=model.rnd_shape(10_000, timegrid.shape[0]))
579 results = model.simulate(timegrid, rnd, expiries=[1.0],
580 power_fwd_prices=[100.0],
581 initial_forecasts={'Region_0': [0.8]})
582 results.get('Power_Germany_FWD0')