Skip to content

Commit

Permalink
Better MC var calculation
Browse files Browse the repository at this point in the history
  • Loading branch information
gavbrennan committed Sep 30, 2022
1 parent 147f515 commit 1da0026
Show file tree
Hide file tree
Showing 5 changed files with 277 additions and 7 deletions.
44 changes: 38 additions & 6 deletions src/Qwack.Models/Risk/McVaRCalculator.cs
Original file line number Diff line number Diff line change
Expand Up @@ -28,49 +28,63 @@ public class McVaRCalculator
private readonly ICurrencyProvider _currencyProvider;
private readonly ICalendarProvider _calendarProvider;
private readonly IFutureSettingsProvider _futureSettingsProvider;
private readonly McModelType _modelType;
private readonly List<IAssetFxModel> _bumpedModels = new();
private readonly Dictionary<string, double> _spotFactors = new();
private readonly Dictionary<string, double[]> _returns = new();
private int _ciIx = 0;

public int CIIX => _ciIx;

public McVaRCalculator(IAssetFxModel model, Portfolio portfolio, ILogger logger, ICurrencyProvider currencyProvider,
ICalendarProvider calendarProvider, IFutureSettingsProvider futureSettingsProvider)
ICalendarProvider calendarProvider, IFutureSettingsProvider futureSettingsProvider, McModelType modelType)
{
_model = model.Clone();
_portfolio = portfolio;
_logger = logger;
_currencyProvider = currencyProvider;
_calendarProvider = calendarProvider;
_futureSettingsProvider = futureSettingsProvider;
_modelType = modelType;
}

public void AddSpotFactor(string assetId, double vol)
{
_spotFactors[assetId] = vol;
}

public void AddReturns(string assetId, double[] returns)
{
_returns[assetId] = returns;
}

public string[] GetSpotFactors() => _spotFactors.Keys.OrderBy(x => x).ToArray();

public void SetCorrelationMatrix(ICorrelationMatrix matrix) => _model.CorrelationMatrix = matrix;

public void CalculateModels()
{
var allAssetIds = _portfolio.AssetIds().Concat(_portfolio.Instruments.Select(x => x.Currency.Ccy).Where(x => x != "USD").Select(x =>$"USD/{x}")).ToArray();
//var allAssetIds = _portfolio.AssetIds().Concat(_portfolio.Instruments.Select(x => x.Currency.Ccy).Where(x => x != "USD").Select(x =>$"USD/{x}")).ToArray();
var allAssetIds = _model.CurveNames.Concat(_portfolio.Instruments.Select(x => x.Currency.Ccy).Where(x => x != "USD").Select(x => $"USD/{x}")).ToArray();
var simulatedIds = allAssetIds.Intersect(_spotFactors.Keys).ToArray();

foreach(var simulatedId in simulatedIds)
{
_model.AddVolSurface(simulatedId, new ConstantVolSurface(_model.BuildDate, _spotFactors[simulatedId]) { AssetId = simulatedId });
if(simulatedId.Length==6 && simulatedId[3]=='/')
_model.FundingModel.VolSurfaces.Add(simulatedId, new ConstantVolSurface(_model.BuildDate, _spotFactors[simulatedId]) { AssetId = simulatedId });
var surf = new ConstantVolSurface(_model.BuildDate, _spotFactors[simulatedId]) { AssetId = simulatedId };
if (_returns.TryGetValue(simulatedId, out var returns))
surf.Returns = returns;
_model.AddVolSurface(simulatedId, surf);
if(simulatedId.Length==6 && simulatedId[3] == '/')
{
_model.FundingModel.VolSurfaces.Add(simulatedId, surf);
}
}

_logger.LogInformation("Simulating {nFac} spot factors", simulatedIds.Length);

var mcSettings = new McSettings
{
McModelType=McModelType.Black,
McModelType= _modelType,
Generator=RandomGeneratorType.MersenneTwister,
NumberOfPaths=2048,
NumberOfTimesteps=2,
Expand Down Expand Up @@ -166,6 +180,24 @@ public Dictionary<string, double> GetContributions(int ix)
return diff.Pivot("TradeId", AggregationAction.Sum).ToDictionary("TradeId").ToDictionary(x => x.Key as string, x => x.Value.Sum(r => r.Value));
}

public decimal ComputeStress(string insId, decimal shockSize)
{
var basePv = _basePvCube.SumOfAllRows;
var baseLevel = _model.GetPriceCurve(insId).GetPriceForFixingDate(_model.BuildDate);
var shockedLevel = baseLevel * Convert.ToDouble(1 + shockSize);

var allScenarios = _resultsCache
.Select(x => (x.Value.SumOfAllRows, _bumpedModels[x.Key].GetPriceCurve(insId).GetPriceForFixingDate(_model.BuildDate)))
.OrderBy(x => x.Item2)
.ToList();

var lr = LinearRegression.LinearRegressionNoVector(allScenarios.Select(x => x.Item2).ToArray(), allScenarios.Select(x => x.SumOfAllRows).ToArray(), false);

var interp = lr.Alpha + lr.Beta * shockedLevel;

return Convert.ToDecimal(interp - basePv);
}

public double CalculateVaR(double ci, Currency ccy) => CalculateVaR(ci, ccy, _portfolio);

private readonly ConcurrentDictionary<int, ICube> _resultsCache = new();
Expand Down
22 changes: 21 additions & 1 deletion src/Qwack.Models/Risk/VaRCalculator.cs
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
using System;
using System.Buffers.Text;
using System.Collections.Concurrent;
using System.Collections.Generic;
using System.Linq;
Expand All @@ -9,6 +10,7 @@
using Qwack.Core.Instruments;
using Qwack.Core.Models;
using Qwack.Math;
using Qwack.Math.Interpolation;
using Qwack.Models.Models;
using Qwack.Models.Risk.Mutators;
using Qwack.Transport.Results;
Expand Down Expand Up @@ -124,7 +126,7 @@ public void AddSurfaceAtmAbsoluteScenarioBumps(string assetId, Dictionary<DateTi

public void CalculateModels()
{
var allAssetIds = _portfolio.AssetIds().Where(x => !(x.Length == 7 && x[3] == '/')).ToArray();
var allAssetIds = _model.CurveNames.Where(x => !(x.Length == 7 && x[3] == '/')).ToArray(); // _portfolio.AssetIds().Where(x => !(x.Length == 7 && x[3] == '/')).ToArray();
var allDatesSet = new HashSet<DateTime>();

if (_spotTypeBumps.Any())
Expand Down Expand Up @@ -290,6 +292,24 @@ public BetaAnalysisResult ComputeBeta(double referenceNav, Dictionary<DateTime,
};
}

public decimal ComputeStress(string insId, decimal shockSize)
{
var basePv = _basePvCube.SumOfAllRows;
var baseLevel = _model.GetPriceCurve(insId).GetPriceForFixingDate(_model.BuildDate);
var shockedLevel = baseLevel * Convert.ToDouble(1 + shockSize);

var allScenarios = _resultsCache
.Select(x => (x.Value.SumOfAllRows, _bumpedModels[x.Key].GetPriceCurve(insId).GetPriceForFixingDate(_model.BuildDate)))
.OrderBy(x=>x.Item2)
.ToList();

var lr = LinearRegression.LinearRegressionNoVector(allScenarios.Select(x => x.Item2).ToArray(), allScenarios.Select(x => x.SumOfAllRows).ToArray(), false);

var interp = lr.Alpha + lr.Beta * shockedLevel;

return Convert.ToDecimal(interp - basePv);
}

public Dictionary<string, double> GetBaseValuations() => _basePvCube.Pivot("TradeId", AggregationAction.Sum).ToDictionary("TradeId").ToDictionary(x => x.Key as string, x => x.Value.Sum(r => r.Value));

public Dictionary<string, double> GetContributions(DateTime scenarioDate)
Expand Down
2 changes: 2 additions & 0 deletions src/Qwack.Options/VolSurfaces/ConstantVolSurface.cs
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,8 @@ public void Build(DateTime originDate, double volatility)
public double InverseCDF(DateTime expiry, double fwd, double p) => VolSurfaceEx.InverseCDFex(this, OriginDate.CalculateYearFraction(expiry, DayCountBasis.Act365F), fwd, p);
public double CDF(DateTime expiry, double fwd, double strike) => this.GenerateCDF2(100, expiry, fwd).Interpolate(strike);

public double[] Returns { get; set; }

public TO_ConstantVolSurface GetTransportObject() => new()
{
AssetId = AssetId,
Expand Down
23 changes: 23 additions & 0 deletions src/Qwack.Options/VolSurfaces/VolSurfaceEx.cs
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,11 @@ public static IInterpolator1D GenerateCDF(this IVolSurface surface, int numSampl

public static IInterpolator1D GenerateCDF2(this IVolSurface surface, int numSamples, DateTime expiry, double fwd, bool returnInverse = false, double strikeScale = 1.0, bool logStrikes = false)
{
if(surface is ConstantVolSurface cv && cv.Returns!=null)
{
return GenerateCDF(cv.Returns, fwd, returnInverse);
}

var premInterp = GeneratePremiumInterpolator(surface, numSamples, expiry, fwd, OptionType.P);
var t = surface.OriginDate.CalculateYearFraction(expiry, DayCountBasis.Act365F);

Expand Down Expand Up @@ -97,6 +102,24 @@ public static IInterpolator1D GenerateCDF2(this IVolSurface surface, int numSamp
InterpolatorFactory.GetInterpolator(x, y, Interpolator1DType.MonotoneCubicSpline);
}

public static IInterpolator1D GenerateCDF(double[] returns, double fwd, bool returnInverse = false)
{
var shockedPrices = returns.Select(r => fwd * (1.0 + r)).OrderBy(x => x).ToArray();
var strikes = shockedPrices.Distinct().ToArray();
var cumProbs = new double[strikes.Length];

var count = 0.0;
for(var i=0;i<strikes.Length;i++)
{
count = shockedPrices.Where(x=>x<=strikes[i]).Count();
cumProbs[i]=count / returns.Length;
}

return returnInverse ?
InterpolatorFactory.GetInterpolator(cumProbs, shockedPrices, Interpolator1DType.MonotoneCubicSpline) :
InterpolatorFactory.GetInterpolator(shockedPrices, cumProbs, Interpolator1DType.MonotoneCubicSpline);
}

public static IInterpolator1D GeneratePDF(this IVolSurface surface, int numSamples, DateTime expiry, double fwd)
{
var deltaK = fwd * 0.0001;
Expand Down
193 changes: 193 additions & 0 deletions src/Qwack.Paths/Processes/TurboSkewSingleAssetFromReturns.cs
Original file line number Diff line number Diff line change
@@ -0,0 +1,193 @@
using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
using Qwack.Core.Models;
using Qwack.Math;
using Qwack.Math.Extensions;
using Qwack.Options.VolSurfaces;
using Qwack.Paths.Features;
using Qwack.Serialization;
using static System.Math;

namespace Qwack.Paths.Processes
{
public class TurboSkewSingleAssetFromReturns : IPathProcess, IRequiresFinish
{
private readonly IATMVolSurface _adjSurface;
private readonly double _correlation;

private readonly DateTime _expiryDate;

private readonly double[] _returns;
private readonly double _atmVol;

private readonly DateTime _startDate;
private readonly int _numberOfSteps;
private readonly string _name;
private readonly Dictionary<DateTime, double> _pastFixings;
private int _factorIndex;
private ITimeStepsFeature _timesteps;
[SkipSerialization]
private readonly Func<double, double> _forwardCurve;
private bool _isComplete;
private double[] _drifts;
private double[] _vols;
private double[] _spotDrifts;
private double[] _spotVols;
private double[] _spotTimesSqrt;
private double _spot0;

private IInterpolator1D[] _invCdfs;

private readonly bool _siegelInvert;

public TurboSkewSingleAssetFromReturns(double[] returns, double atmVol, DateTime startDate, DateTime expiryDate, int nTimeSteps, Func<double, double> forwardCurve, string name, Dictionary<DateTime, double> pastFixings = null, IATMVolSurface fxAdjustSurface = null, double fxAssetCorrelation = 0.0)
{
_returns = returns;
_atmVol = atmVol;

_startDate = startDate;
_expiryDate = expiryDate;
_numberOfSteps = nTimeSteps;
_name = name;
_forwardCurve = forwardCurve;
_pastFixings = pastFixings ?? (new Dictionary<DateTime, double>());

_adjSurface = fxAdjustSurface;
_correlation = fxAssetCorrelation;
}

public bool IsComplete => _isComplete;

public void Finish(IFeatureCollection collection)
{
if (!_timesteps.IsComplete)
{
return;
}

//drifts and vols...
_drifts = new double[_timesteps.TimeStepCount];
_vols = new double[_timesteps.TimeStepCount];
_spotDrifts = new double[_timesteps.TimeStepCount];
_spotVols = new double[_timesteps.TimeStepCount];
_invCdfs = new IInterpolator1D[_timesteps.TimeStepCount];

_spot0 = _forwardCurve(0);
var prevSpot = _spot0;
for (var t = 1; t < _drifts.Length; t++)
{
var atmVol = _atmVol;
var fxAtmVol = _adjSurface == null ? 0.0 : _adjSurface.GetForwardATMVol(0, _timesteps.Times[t]);
var driftAdj = _adjSurface == null ? 1.0 : Exp(atmVol * fxAtmVol * _timesteps.Times[t] * _correlation);
var spot = _forwardCurve(_timesteps.Times[t]) * driftAdj;
var varStart = Pow(atmVol, 2) * _timesteps.Times[t - 1];
var varEnd = Pow(atmVol, 2) * _timesteps.Times[t];
var fwdVariance = Max(0, varEnd - varStart);
_vols[t] = Sqrt(fwdVariance / _timesteps.TimeSteps[t]);
_drifts[t] = Log(spot / prevSpot) / _timesteps.TimeSteps[t];
_invCdfs[t] = VolSurfaceEx.GenerateCDF(_returns, spot, true);

_spotVols[t] = atmVol;
_spotDrifts[t] = Log(spot / _spot0) / _timesteps.Times[t];

prevSpot = spot;
}

_spotTimesSqrt = _timesteps.Times.Select(x => Sqrt(x)).ToArray();

_isComplete = true;
}

public void Process(IPathBlock block)
{
for (var path = 0; path < block.NumberOfPaths; path += Vector<double>.Count)
{
var previousStep = new Vector<double>(_forwardCurve(0));
var steps = block.GetStepsForFactor(path, _factorIndex);
var c = 0;
foreach (var kv in _pastFixings.Where(x => x.Key < _startDate))
{
steps[c] = new Vector<double>(kv.Value);
c++;
}
steps[c] = previousStep;

if (_siegelInvert)
{
for (var step = c + 1; step < block.NumberOfSteps; step++)
{
var W = steps[step];
var dt = new Vector<double>(_timesteps.TimeSteps[step]);
var bm = (_drifts[step] + _vols[step] * _vols[step] / 2.0) * dt + (_vols[step] * _timesteps.TimeStepsSqrt[step] * W);
previousStep *= bm.Exp();
steps[step] = previousStep;
}

//transform
for (var step = c + 1; step < block.NumberOfSteps; step++)
{
var ws = new double[Vector<double>.Count];
for (var v = 0; v < ws.Length; v++)
{
var t1 = Log(steps[step][v] / _spot0);
var t2 = (_spotDrifts[step] + _spotVols[step] * _spotVols[step] / 2.0) * _timesteps.Times[step];
t1 -= t2;
t1 /= _spotVols[step] * _spotTimesSqrt[step];
t1 = Statistics.NormSDist(t1);
ws[v] = _invCdfs[step].Interpolate(t1);
}
steps[step] = new Vector<double>(ws);
}
}
else
{
for (var step = c + 1; step < block.NumberOfSteps; step++)
{
var W = steps[step];
var dt = new Vector<double>(_timesteps.TimeSteps[step]);
var bm = (_drifts[step] - _vols[step] * _vols[step] / 2.0) * dt + (_vols[step] * _timesteps.TimeStepsSqrt[step] * W);
previousStep *= bm.Exp();
steps[step] = previousStep;
}

//transform
for (var step = c + 1; step < block.NumberOfSteps; step++)
{
var ws = new double[Vector<double>.Count];
for (var v = 0; v < ws.Length; v++)
{
var t1 = Log(steps[step][v] / _spot0);
var t2 = (_spotDrifts[step] - _spotVols[step] * _spotVols[step] / 2.0) * _timesteps.Times[step];
t1 -= t2;
t1 /= _spotVols[step] * _spotTimesSqrt[step];
t1 = Statistics.NormSDist(t1);
ws[v] = _invCdfs[step].Interpolate(t1);
}
steps[step] = new Vector<double>(ws);
}
}
}
}

public void SetupFeatures(IFeatureCollection pathProcessFeaturesCollection)
{
var mappingFeature = pathProcessFeaturesCollection.GetFeature<IPathMappingFeature>();
_factorIndex = mappingFeature.AddDimension(_name);

_timesteps = pathProcessFeaturesCollection.GetFeature<ITimeStepsFeature>();
_timesteps.AddDates(_pastFixings.Keys.Where(x => x < _startDate));

var stepSize = (_expiryDate - _startDate).TotalDays / _numberOfSteps;
var simDates = new List<DateTime>();
for (var i = 0; i < _numberOfSteps; i++)
{
simDates.Add(_startDate.AddDays(i * stepSize).Date);
}

_timesteps.AddDates(simDates.Distinct());
}
}
}

0 comments on commit 1da0026

Please sign in to comment.