Correlação entre rendibilidade de activos

math = require('mathjs')
//import {multivariateNormal} from '@sw1227/multivariate-normal-distribution'
import {Plot} from "@mkfreeman/plot-tooltip"
import {Scrubber} from "@mbostock/scrubber"
import {Legend, Swatches} from "@d3/color-legend"

function choleskyDecomposition(matrix) {
  // Argument "matrix" can be either math.matrix or standard 2D array
  const A = math.matrix(matrix);
  // Matrix A must be symmetric
  console.assert(math.deepEqual(A, math.transpose(A)));

  const n = A.size()[0];
  // Prepare 2D array with 0
  const L = new Array(n).fill(0).map(_ => new Array(n).fill(0));

  d3.range(n).forEach(i => {
    d3.range(i+1).forEach(k => {
      const s = d3.sum(d3.range(k).map(j => L[i][j]*L[k][j]));
      L[i][k] = i === k ? math.sqrt(A.get([k, k]) - s) : 1/L[k][k] * (A.get([i, k]) - s);
    })
  });
  return L;
}

// Generate n samples from standard normal distribution
// [Box-Muller transform / sw1227 | Observable](https://observablehq.com/@sw1227/box-muller-transform)
function boxMuller(n) {
  const samples = [];
  Array(Math.ceil(n / 2)).fill().forEach(_ => {
    const R = Math.sqrt(-2 * Math.log(Math.random()));
    const theta = 2 * Math.PI * Math.random();
    samples.push(R * Math.cos(theta)); // z1
    samples.push(R * Math.sin(theta)); // z2
  });
  // if n is odd, drop the last element
  return samples.slice(0, n);
}

// [Multivariate Normal Distribution / sw1227 | Observable](https://observablehq.com/@sw1227/multivariate-normal-distribution)
function multivariateNormal(mean, covArray) {
  const n = mean.length;
  const cov = math.matrix(covArray);
  return {
    // Probability Density Function
    pdf: x => {
      const c = 1 / (math.sqrt(2*math.PI)**n * math.sqrt(math.det(cov)));
      return c * math.exp(
        -(1/2) * math.multiply(
          math.subtract(math.matrix(x), math.matrix(mean)),
          math.inv(cov),
          math.subtract(math.matrix(x), math.matrix(mean))
        )
      );
    },
    // Differential entropy
    entropy: 0.5*math.log(math.det(cov)) + 0.5*n*(1 + math.log(2*math.PI)),
    // Generate n samples using Cholesky Decomposition
    sample: n_samples => Array(n_samples).fill().map(_ => {
      const L = choleskyDecomposition(cov);
      const z = boxMuller(n);
      return math.add(
        math.matrix(mean),
        math.multiply(L, math.matrix(z))
      ).toArray();
    }),
  };
}

Considere dois activos X e Y com rendibilidades que seguem uma Distribuição Normal.

As rendibilidades esperadas, desvios-padrão e correlação entre as rendibilidades dos activos são dadas por:

viewof x_mean = Inputs.range([0, 0.4], {value: 0.01, step: 0.01, label: "Rendibilidade esperada de X", width: "150px"})
viewof y_mean = Inputs.range([0, 0.4], {value: 0.02, step: 0.01, label: "Rendibilidade esperada de Y", width: "150px"})
viewof x_sd = Inputs.range([0.01, 1], {value: 0.15, step: 0.01, label: "Desvio padrão das rendibilidades de X", width: "150px"})
viewof y_sd = Inputs.range([0.01, 1], {value: 0.20, step: 0.01, label: "Desvio padrão das rendibilidades de Y", width: "150px"})
viewof corr = Inputs.range([-1, 1], {value: 0.0, step: 0.01, label: "Correlação entre as rendibilidades de X e de Y", width: "150px"})

Nota: Sempre que alterar algum dos parâmetros uma nova amostra é gerada.

cov = math.matrix([[Math.pow(x_sd,2), corr*x_sd*y_sd], [corr*x_sd*y_sd, Math.pow(y_sd,2)]])
// cov

// Generate a multivariate normal random variates
normal = multivariateNormal([x_mean, y_mean], cov)
// Extract the random variates
random_variates = normal.sample(50) 
sample_returns = random_variates.map(d => ({x: d[0], y: d[1]}))

sample_prices = {
  var sample_prices = [];
  var sums = [0, 0];
  var count = 0;
  for(var i in sample_returns){
    count += 1;
    sums[0] += sample_returns[i].x;
    sums[1] += sample_returns[i].y;
    sample_prices.push({
      i: count, 
      price_x: 100*Math.exp(sums[0]),
      price_y: 100*Math.exp(sums[1])
    });
  }
  return sample_prices;
}
//sample_prices
//sample_prices = sample_returns.map(d => {prices_x: cumulativeSum(x)})

//sample_returns
//Inputs.table(sample_returns)

iterator = Array.from({length: sample_returns.length}, (_, i) => i+1)
viewof i = Scrubber(iterator, {delay: 700})
Plot.plot({
  marks: [
    Plot.dot(sample_returns.slice(0, i), {
      x: "x", 
      y: "y", 
      fill: "orange",
      title: (d) =>
        `Rendibilidades \n X: ${d.x.toFixed(3)} \n Y: ${d.y.toFixed(3)}`
    }),
    Plot.dot(sample_returns.slice(i-1, i), {
      x: "x", 
      y: "y", 
      fill: "green",
      title: (d) =>
        `Rendibilidades \n X: ${d.x.toFixed(3)} \n Y: ${d.y.toFixed(3)}`
    })
  ]
})
Plot.plot({
  marks: [
    Plot.line(sample_prices.slice(0, i), {x: "i", y: "price_x", stroke: "blue"}),
    Plot.dot(sample_prices.slice(0, i), {
      x: "i", 
      y: "price_x", 
      stroke: "blue",
      title: (d) =>
        `Tempo ${d.i.toFixed(0)} \n Preço de X: ${d.price_x.toFixed(3)}`
    }),
    Plot.dot(sample_prices.slice(i-1, i), {
      x: "i", 
      y: "price_x", 
      fill: "green",
      title: (d) =>
        `Tempo ${d.i.toFixed(0)} \n Preço de X: ${d.price_x.toFixed(3)}`
    }),
    Plot.line(sample_prices.slice(0, i), {x: "i", y: "price_y", stroke: "red"}),
    Plot.dot(sample_prices.slice(0, i), {
      x: "i", 
      y: "price_y", 
      stroke: "red",
      title: (d) =>
        `Tempo ${d.i.toFixed(0)} \n Preço de Y: ${d.price_y.toFixed(3)}`
    }),
    Plot.dot(sample_prices.slice(i-1, i), {
      x: "i", 
      y: "price_y", 
      fill: "green",
      title: (d) =>
        `Tempo ${d.i.toFixed(0)} \n Preço de Y: ${d.price_y.toFixed(3)}`
    }),
    Plot.ruleX([0]), // a rule at x = 0
    Plot.ruleY([0])
  ]
})
Swatches(d3.scaleOrdinal(["X", "Y"], ["blue", "red"]))

O gráfico abaixo mostra 50 rendibilidades simuladas destes dois activos:

O gráfico abaixo mostra os preços simulados desses dois activos, tratando as rendibilidades simuladas como rendibilidades logarítmicas (em capitalização contínua) e assumindo que os seus preços no momento zero são de 100: