fix: portfolio volatility calculation

allocation/risk.go

@@ -36,7 +36,7 @@ func (*EqualRiskContribution) PolicyWeights(ctx context.Context, _ time.Time, as
 	copy(weights, ws)
 
 	return weights, optWeights(ctx, weights, func(ws []float64) float64 {
-		_, _, riskWeights := calculations.RiskFromRiskContribution(assetRisks, ws, cm)
+		riskWeights := calculations.RiskWeights(calculations.PortfolioVolatility(assetRisks, ws, cm))
 		var diff float64
 		for i := range riskWeights {
 			diff += math.Abs(target - riskWeights[i])

calculations/correlation.go

@@ -0,0 +1,29 @@
+package calculations
+
+import (
+	"gonum.org/v1/gonum/stat"
+)
+
+func CorrelationMatrix(values [][]float64) [][]float64 {
+	n := len(values)
+	if len(values) == 0 {
+		return nil
+	}
+	m := make([][]float64, n)
+	for i := range m {
+		m[i] = make([]float64, n)
+	}
+	for i := 0; i < n; i++ {
+		for j := i; j < n; j++ {
+			if i == j {
+				m[i][j] = 1.0
+				continue
+			}
+			vi, vj := values[i], values[j]
+			c := stat.Correlation(vi, vj, nil)
+			m[i][j] = c
+			m[j][i] = c
+		}
+	}
+	return m
+}

calculations/risk.go

@@ -2,117 +2,12 @@ package calculations
 
 import (
 	"errors"
-	"fmt"
 	"math"
+	"slices"
 
-	"gonum.org/v1/gonum/floats"
-	"gonum.org/v1/gonum/mat"
 	"gonum.org/v1/gonum/stat"
 )
 
-func DenseMatrixToFloatSlices(dense *mat.Dense) [][]float64 {
-	if dense == nil {
-		return make([][]float64, 0)
-	}
-	iLength, jLength := dense.Dims()
-	result := make([][]float64, iLength)
-	d := dense.RawMatrix().Data
-	for i := range result {
-		result[i] = d[i*jLength : (i+1)*jLength]
-	}
-	return result
-}
-
-func Risk(risks, weights []float64, correlations [][]float64) (float64, error) {
-	numberOfAssets := len(risks)
-
-	for i, r := range risks {
-		risks[i] = r / 100.0
-	}
-	if sum := floats.Sum(weights); math.Abs(100.0-sum) > 0.0001 {
-		return 0, fmt.Errorf(`weights must add up to 100 but got %.2f`, sum)
-	}
-	for i, w := range weights {
-		weights[i] = w / 100.0
-	}
-
-	if len(weights) != numberOfAssets {
-		return 0, errors.New("the number of risks must be the same as the number of weights")
-	}
-
-	if len(correlations) != numberOfAssets {
-		return 0, fmt.Errorf("correlations must be n by n where n is the length of risks; the number of rows is should be %d", numberOfAssets)
-	}
-
-	cor := mat.NewDense(len(risks), len(risks), nil)
-	for i, row := range correlations {
-		if len(row) != numberOfAssets {
-			return 0, fmt.Errorf("correlations must be n by n where n is the length of risks; row %d has %d it should have %d", i, len(row), numberOfAssets)
-		}
-		for j, v := range correlations[i] {
-			if v < -1 || v > 1 {
-				return 0, fmt.Errorf("correlation values must be within [-1, 1] the value %.2f at row %d column %d is not in this range", v, i+1, j+1)
-			}
-			cor.Set(i, j, v)
-		}
-	}
-
-	r, _, _ := RiskFromRiskContribution(risks, weights, cor)
-
-	return r, nil
-}
-
-func RiskFromRiskContribution(risks, weights []float64, correlations *mat.Dense) (float64, []float64, []float64) {
-	n := len(risks)
-	nSquared := n * n
-
-	// START memory allocation stuff
-	bufferSize := 3*nSquared + 3*n
-	b := make([]float64, bufferSize)
-	b1 := b[:nSquared:nSquared]
-	b = b[len(b1):]
-	b2 := b[:nSquared:nSquared]
-	b = b[len(b2):]
-	b3 := b[:n:n]
-	b = b[len(b3):]
-	b4 := b[:nSquared:nSquared]
-	b = b[len(b4):]
-	v := b[:n:n]
-	b = b[len(v):]
-	rw := b[:n:n]
-	// END memory allocation stuff
-
-	d := mat.NewDense(n, n, b1)
-
-	for i := 0; i < len(risks); i++ {
-		d.Set(i, i, risks[i])
-	}
-
-	V := mat.NewDense(n, n, b2)
-	V.Product(d, correlations, d)
-
-	mx1 := mat.NewDense(n, 1, b3)
-	mx1.Product(V, mat.NewDense(n, 1, weights))
-	mx2 := mat.NewDense(n, n, b4)
-	mx2.Product(mx1, mat.NewDense(1, n, weights))
-
-	for i := 0; i < n; i++ {
-		v[i] = mx2.At(i, i)
-	}
-	pv := floats.Sum(v)
-	copy(rw, v)
-	for i := range rw {
-		rw[i] /= pv
-	}
-
-	return math.Sqrt(pv), v, rw
-}
-
-func ExpectedRisk(risks, weights []float64, correlations *mat.Dense) float64 {
-	r, _, _ := RiskFromRiskContribution(risks, weights, correlations)
-	return r
-}
-
 func NumberOfBets(weightedAverageRisk float64, portfolioRisk float64) (float64, error) {
 	if portfolioRisk == 0 {
 		return 0, errors.New("can't divide by 0")
@@ -136,3 +31,50 @@ func RiskFromStdDev(values []float64) float64 {
 func WeightedAverageRisk(weights, risks []float64) float64 {
 	return stat.Mean(risks, weights)
 }
+
+func RiskWeights(portfolioVol float64, riskContributions []float64) []float64 {
+	result := slices.Clone(riskContributions)
+	for i := range result {
+		result[i] = result[i] / portfolioVol
+	}
+	return result
+}
+
+func PortfolioVolatility(weights, stdDevs []float64, correlations [][]float64) (float64, []float64) {
+	var variance float64
+	n := len(weights)
+	covariances := calculateCovarianceMatrix(weights, stdDevs, correlations)
+
+	for i := 0; i < n; i++ {
+		for j := 0; j < n; j++ {
+			variance += covariances[i][j]
+		}
+	}
+
+	totalRisk := math.Sqrt(variance)
+	riskContributions := make([]float64, n)
+
+	for i := 0; i < n; i++ {
+		for j := 0; j < n; j++ {
+			riskContributions[i] += covariances[i][j]
+		}
+		riskContributions[i] = riskContributions[i] / totalRisk
+	}
+
+	return totalRisk, riskContributions
+}
+
+func calculateCovarianceMatrix(weights, vols []float64, correlations [][]float64) [][]float64 {
+	n := len(weights)
+	covariances := make([][]float64, n)
+	for i := range covariances {
+		covariances[i] = make([]float64, n)
+	}
+
+	for i := 0; i < n; i++ {
+		for j := 0; j < n; j++ {
+			covariances[i][j] = weights[i] * weights[j] * vols[i] * vols[j] * correlations[i][j]
+		}
+	}
+	return covariances
+}

calculations/risk_test.go

@@ -2,34 +2,14 @@ package calculations
 
 import (
 	"fmt"
-	"math"
 	"testing"
 
 	"github.com/portfoliotree/round"
 	"github.com/stretchr/testify/assert"
-
-	"gonum.org/v1/gonum/floats"
-	"gonum.org/v1/gonum/mat"
 )
 
 func TestRisk(t *testing.T) {
-	t.Run("simple", func(t *testing.T) {
-		risk, err := Risk(
-			[]float64{10, 10},
-			[]float64{100, 0.0},
-			[][]float64{
-				{1, 1},
-				{1, 1},
-			},
-		)
-		if err != nil {
-			t.Errorf("it should not return an error: %s", err)
-		}
-
-		if math.IsNaN(risk) {
-			t.Errorf("risk should be a number got: %f", risk)
-		}
-	})
+
 }
 
 func ExampleRiskFromStdDev() {
@@ -108,26 +88,63 @@ func TestNumberOfBets_Zero_Portfolio_Risk(t *testing.T) {
 	}
 }
 
-func TestRiskFromRiskContribution(t *testing.T) {
-	ws := []float64{-.1, .8, .3}
-	rs := []float64{.2, .05, .25}
-	cs := mat.NewDense(3, 3, []float64{
-		1.00, 0.15, 0.60,
-		0.15, 1.00, 0.05,
-		0.60, 0.05, 1.00,
-	})
-
-	const precision = 4
-	tr, rcs, rw := RiskFromRiskContribution(rs, ws, cs)
-	tr = round.Decimal(tr, precision)
-	_ = round.Recursive(rcs, precision)
-	_ = round.Recursive(rw, precision)
-	sum := round.Decimal(floats.Sum(rw), precision)
-
-	assert.Equal(t, rcs, []float64{-0.0006, 0.0016, 0.0049})
-	assert.Equal(t, tr, .0767)
-	assert.Equal(t, rw, []float64{-0.1054, 0.2770, 0.8284})
-	assert.Equal(t, sum, 1.0)
+func TestPortfolioVolatility(t *testing.T) {
+	for _, tt := range []struct {
+		Name         string
+		Weights      []float64
+		Risks        []float64
+		Correlations [][]float64
+
+		ExpectedPortfolioVolatility float64
+		ExpectRiskContributions     []float64
+	}{
+		{
+			Name:    "legacy example",
+			Weights: []float64{-.1, .8, .3},
+			Risks:   []float64{.2, .05, .25},
+			Correlations: [][]float64{
+				{1.00, 0.15, 0.60},
+				{0.15, 1.00, 0.05},
+				{0.60, 0.05, 1.00},
+			},
+			ExpectRiskContributions:     []float64{-0.0081, 0.0212, 0.0635},
+			ExpectedPortfolioVolatility: .0767,
+		},
+		{
+			Name:    "two",
+			Weights: []float64{0.5, 0.5},
+			Risks:   []float64{0.1668, 0.0428},
+			Correlations: [][]float64{
+				{1.00, 0.25},
+				{0.25, 1.00},
+			},
+			ExpectRiskContributions:     []float64{0.0812, 0.0099},
+			ExpectedPortfolioVolatility: .0911,
+		},
+		{
+			Name:    "three",
+			Weights: []float64{0.2500, 0.5000, 0.2500},
+			Risks:   []float64{0.0428, 0.1668, 0.1693},
+			Correlations: [][]float64{
+				{1.0000, 0.2465, 0.4077},
+				{0.2465, 1.0000, 0.1091},
+				{0.4077, 0.1091, 1.0000},
+			},
+			ExpectRiskContributions:     []float64{0.0051, 0.0740, 0.0231},
+			ExpectedPortfolioVolatility: .1022,
+		},
+	} {
+		t.Run(tt.Name, func(t *testing.T) {
+			portfolioVol, riskContributions := PortfolioVolatility(tt.Weights, tt.Risks, tt.Correlations)
+
+			const precision = 4
+			portfolioVol = round.Decimal(portfolioVol, precision)
+			_ = round.Recursive(riskContributions, precision)
+
+			assert.Equal(t, tt.ExpectRiskContributions, riskContributions)
+			assert.Equal(t, tt.ExpectedPortfolioVolatility, portfolioVol)
+		})
+	}
 }
 
 // func TestLegacyRiskFromRiskContribution(t *testing.T) {