• Amplamente utilizadas nos mercados financeiros
• A implementação dessa estratégia de arbitragem estatística pode restringir a especulação excessiva e aumentar a liquidez do mercado.
• Essa técnica consiste em uma operação onde o investidor combina uma posição comprada (long position) com uma venda a descoberto (short position) de dois ativos extremamente correlacionados.
• Quando o equilíbrio é novamente restaurado, ambas as posições são desfeitas e o investidor obtém lucro.

• Movimento Browniano
• Movimento Browniano geométrico
• Modelo Black and Scholes
• Ornstein–Uhlenbeck
• Cox-Ingersoll-Ross

Também chamado de processo de Wiener, \(W = \{W(t), t\geq0\}\), é um processo Gaussiano com caminhos contínuos e com incrementos independentes tais que \(W(0) = 0\) com probabilidade \(1\), \(\mathbb{E}W(t) = 0\) e \(\mathrm{Var}(W(t) - W(s)) = t - s\), para todo \(0\leq s \leq t\). Na prática, o que é relevante para nossos propósitos é que \(W(t) − W(s) \sim N(0, t − s)\), para \(0\leq s \leq t\) e que em quaisquer dois intervalos disjuntos, digamos \((t_1, t_2)\), \((t_3, t_4)\) com \(t1 \leq t2 \leq t3 \leq t4\), o os incrementos \(W(t_2) − W(t_1)\) e \(W(t_4) − W(t_3)\) são independentes.

Dado um incremento de tempo fixo \(\Delta t >0\), pode-se facilmente simular uma trajetória do processo de Wiener no intervalo de tempo \([0, T]\). De fato, para \(W_{\Delta t}\) é verdade que: \[ W(\Delta t) = W(\Delta t) - W(0) \sim N(0, \Delta t) \sim \sqrt{\Delta t} N (0,1), \]

e o mesmo vale para qualquer outro incremento \(W(t + \Delta t) - W(t)\): \[ W(t+ \Delta t) - W(t) \sim N(0, \Delta t) \sim \sqrt{\Delta t} N (0,1), \]

Assim, podemos simular um desses caminhos da seguinte maneira. Divida o intervalo \([0, T]\) em uma grade como \(0 = t_1 < t_2 < \dots < t_{N−1} < t_N = T\) com \(t_{i+1} − t_i = \Delta t\). Defina \(i = 1\) e \(W(0) = W(t_1) = 0\) e itere o seguinte algoritmo.

• Gera um (novo) número aleatório \(z\) a partir da distribuição gaussiana padrão.
• \(i = i + 1\).
• Set \(W(t_i) = W(t_{i-1}) + z\sqrt{\Delta t}\)
• Se \(i \leq N\), itera a partir do passo 1

Traduzindo este algoritmo acima como codificação.

set.seed(123)
N <- 100 
T <- 1 
Delta <- T/ N 
W <- numeric ( N +1) 
t <- seq (0 ,T , length = N +1)
for (i in 2:( N +1)){
  W [i ] <- W [i -1] + rnorm (1) * sqrt ( Delta )
}
plot(t ,W , type = "l " , main = "Wiener process" , ylim = c ( -1 ,1))

Este processo tem a propriedade de ter incrementos multiplicativos independentes e é definido como uma função do movimento Browniano padrão \[ S(t) = x \mathrm{\exp}\left\{\left(r - \frac{\sigma^2}{2}\right)t + \sigma W(t) \right\}, \quad t>0, \]

com \(S(0) = x, \quad x \in \mathbb{R}\) é o valor inicial, \(\sigma>0\) interpretado como a volatilidade e \(r\) interpretado como a taxa de juros. Este processo foi introduzido no contexto financeiro em Osborne (1959). Uma implementação para simular um caminho do movimento Browniano geométrico pode ser baseada em algoritmos anteriores. Assumindo que \(W\) contém uma trajetória de \(W\) e \(t\) é o vetor contendo todos os pontos no tempo, um caminho pode ser traçado da seguinte forma.

Uma forma equivalente de simular uma trajetória do movimento Browniano geométrico é simulando os incrementos de \(S\). \[ \scriptstyle S(t+\Delta t) = S(t) \mathrm{\exp}\left\{\left(r - \frac{\sigma^2}{2}\right)\Delta t + \sigma\sqrt{ \Delta t}Z\right\}, \] com \(Z \sim N(0,1)\). A equação é um caso particular do movimento Browniano geométrico generalizado, que é um processo partindo de \(x\) no tempo \(s\) cuja dinâmica é \[ \scriptstyle Z_{s,x}(t) = x \mathrm{\exp}\left\{\left(r - \frac{\sigma^2}{2}\right)(t-s) + \sigma (W(t ) - W(s))\right\}, \quad t>s, \]

Uma ferramenta importante do cálculo estocástico que também é útil em simulações é a fórmula de Itô. Esta fórmula pode ser vista como a versão estocástica de uma expansão de Taylor de \(f(X)\) parada na segunda ordem, onde \(X\) é um processo de difusão. O lema de Itô diz que se \(f(t, x)\) é uma função duas vezes diferenciável tanto em \(t\) quanto em \(x\), então \[ \scriptstyle f(t, X_t) = f(0, X_0) + \int_{0}^{t}f_t(u, X_u)du + \int_{0}^{t}f_x(u, X_u)dX_u + \dfrac{1}{2}\int_{0}^{t}f_{xx}(u, X_u)(dX_u)^2, \]

em que

\[ \scriptstyle \dfrac{\partial f(t,x)}{\partial t} = f_t(t, x), \quad \dfrac{\partial f(t,x)}{\partial x} = f_x(t, x), \quad \dfrac{\partial^2 f(t,x)}{\partial x^2} = f_{xx}(t, x) \] ou, na forma diferencial \[ \scriptstyle df(t, X_t) = f_t(t, X_t)dt + f_x(t, X_t)dX_t + \dfrac{1}{2}f_{xx}(t, X_t)(dX_t)^2. \]

Após as duas definições acima de Black e Scholes, ou o modelo de movimento Browniano geométrico resolve a equação diferencial estocástica \[ \scriptstyle dX_t = \theta_1X_tdt + \theta_2X_tdW_t, \quad X_0 = x_0, \quad \theta_1 \in \mathbb{R},\theta_2 \in \mathbb{R}_+. \] com \(\theta_2>0\). O parâmetro \(\theta_1\) é interpretado como a taxa de juros constante e \(\theta_2\) como a volatilidade das atividades de risco. A solução explícita é

\[ \scriptstyle X_t = x_0e^{(\theta_1 - \frac{1}{2}\theta^2_2)t + \theta_2W_t} \]

A função de densidade condicional é log-normal com a média e a variância de sua transformada logarítmica (isto é, a média logarítmica e a variância logarítmica) dada por \[ \scriptstyle \mu = log(x_0) + \left(\theta_1 - \frac{1}{2}\theta^2_2\right)t, \qquad \sigma^2 = \theta^2_2t \]

com média e variância

Por isso

\[\scriptstyle p_{\theta}(t, y|x_0) = \dfrac{1}{\theta_{2}y\sqrt{2\pi t}}\exp{\left\{-\dfrac{(\log{y} - (\log {x_0} + (\theta_1 -\frac{1}{2}\theta^2_2)t))^2}{2\theta^{2}_2t}\right\}} \]

O processo Ornstein-Uhlenbeck ou Vasicek é a única solução para a seguinte equação diferencial estocástica \[ \scriptstyle dX_t = (\theta_1 - \theta_2X_t)dt + \theta_3dW_t, \quad X_0 = x_0, \] com \(\theta_1, \theta_2 \in \mathbb{R}, \theta3 \in \mathbb{R}+\) e é ergódico para \(\theta_2 > 0\). Para \(\theta_2 > 0\), este é um processo de reversão à média, o que significa que o processo tende a oscilar em torno de algum estado de equilíbrio. Outra propriedade interessante do processo de Ornstein-Uhlenbeck é que, ao contrário do movimento Browniano, é um processo com variância finita para tudo \(t \geq 0\). Outra parametrização do processo Ornstein-Uhlenbeck mais comum em modelagem financeira é \[ dX_t = \theta(\mu - X_t)dt + \sigma dW_t, \quad X_0 = x_0, \]

onde \(\sigma\) é interpretado como a volatilidade, \(\mu\) é o valor de equilíbrio de longo prazo do processo e \(\theta\) é a velocidade de reversão. Como aplicação do lema de Itô, podemos mostrar a solução explícita de escolhendo \(f(t, x) = xe^{\theta_2 t}\).

\[ \scriptstyle X_t = \dfrac{\theta_1}{\theta_2} + \left(x_0 - \dfrac{\theta_1}{\theta_2}\right)e^{-\theta_2 t} + \theta_3\int_{0}^{t}e^{-\theta_2(t-u)}dW_u \] ou, na segunda parametrização, \[ \scriptstyle X_t = \mu + (x_0 - \mu)e^{-\theta t} + \sigma\int_{0}^{t}e^{-\theta(t-u)}dW_u \]

A lei invariante

Para \(\theta_2 > 0\), o processo também é ergódico, e sua lei invariante é a densidade gaussiana com média \(\theta_1/\theta_2\) e variância \(\theta^{2}_3/2\theta_2\), \[ \scriptstyle X_t \sim N\left(\frac{\theta_1}{\theta_2}, \frac{\theta^{2}_3}{2\theta_2}\right). \]

A lei condicional

Para qualquer \(t \geq 0\), o processo de Ornstein-Uhlenbeck tem uma densidade de transição Gaussiana (ou condicional) \(p_{\theta}(t, X_t|X_0=x_0)\); ou seja, a densidade da distribuição de \(X_t\) dado \(X_0 = x_0\), com média e variância respectivamente

\[ \scriptstyle m(t,x) = \mathbb{E}_\theta(X_t | X_0 = x_0) = \dfrac{\theta_1}{\theta_2} + \left(x_0 - \frac{\theta_1}{\theta_2}\right)e^{-\theta_2 t} \] e \[ \scriptstyle \nu(t,x) = \mathrm{Var}_\theta(X_t | X_0 = x_0) = \dfrac{\theta^{2}_3(1 - e^{-2\theta_2 t})}{2\theta_2}. \]

Exemplo de simulação da integral estocástica

Caminho simulado do processo Ornstein-Uhlenbeck \(dX_t = −\theta X_t + \sigma dWt\) com \(X(0) = 10\), \(\theta = 5\) e \(\sigma = 3,5\)

library(sde)
W<-BM()
t<-time(W)
N<-length(t)
x<-10
theta<-5
sigma<-3.5
X<-numeric(N)
X [1]<-x
ito.sum<-c(0, sapply (2: N , function(x){
  exp(-theta*(t[x] -t[x-1]))*(W[x]- W[x-1])
  }))
X<-sapply(1:N , function(x){X[1]*exp(-theta*t[x]) + sum(ito.sum[1:x])})
X<-ts(X , start = start(W), deltat = deltat(W))
plot(X, main = "Processo Ornstein - Uhlenbeck")

Uma outra família interessante de modelos paramétricos é a do processo Cox-Ingersoll-Ross. Este modelo foi introduzido por Feller como um modelo para crescimento populacional (ver Feller (1951) e Feller et al. (1951)) e tornou-se bastante popular em finanças depois que Cox, Ingersoll e Ross o propuseram para modelar taxas de juros de curto prazo Cox, Ingersoll, and Ross (1985). O processo CIR é a solução para a equação diferencial estocástica.

\[ \scriptstyle dX_t = (\theta_1 - \theta_2X_t)dt + \theta_3\sqrt{X_t}dW_t, \quad X_0 = x_0 > 0, \]

às vezes parametrizado como

\[ \scriptstyle dX_t = \theta(\beta - \theta_2X_t)dt + \sigma\sqrt{X_t}dW_t, \quad X_0 = x_0 > 0, \] onde \(\theta_1, \theta_2, \theta_3 \in \mathbb{R}_+\). Se \(2\theta_1 > \theta^2_3\), o processo é estritamente positivo caso contrário é não negativo, o que significa que pode atingir o estado \(0\). A equação diferencial estocástica tem a solução explícita

\[ \scriptstyle X_t = \left(X_0 - \dfrac{\theta_1}{\theta_2}\right)e^{-\theta_2 t} + \theta_3 e^{-\theta_2 t}\int_{0}^{t}e^{-\theta_2 t}\sqrt{X_u}dW_u \]

• Método de máxima verossimilhança.
• Em Pyke (1963), Rao (1972), Rao (1983) e Dacunha-Castelle and Florens-Zmirou (1986), são dadas condições suficientes para a consistência e normalidade assintótica desses estimadores de máxima verossimilhança.
• Em particular, em Dacunha-Castelle and Florens-Zmirou (1986) a consistência e a normalidade assintótica são provadas, independentemente do tamanho de \(\Delta\), o que parece ser uma boa propriedade em aplicações.
• Função mle() do pacote stats4
• Também foi usado funções do pacote sde para a simulação dos modelos apresentados.

library(sde)
library(forecast)
#LogLik BS
BS.lik <- function(theta, sigma){
  n <- length(X)
  dt <- deltat(X)
  -sum(sde::dcBS(x = X[2:n], Dt = dt , x0 = X[1:(n -1)], 
                 theta = c(theta, sigma), log = TRUE))
}
set.seed(123)
X<-sde.sim(model = "BS", theta = c(0.5, 0.2), N = 1000, delta = 0.01)
autoplot(X, main = 'Simulação para o modelo Black e Scholes')
fit<-mle(BS.lik, start = list(theta =1 , sigma =1),
         method = "L-BFGS-B" , lower = c(0.01 ,0.01)) 
summary(fit)

## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = BS.lik, start = list(theta = 1, sigma = 1), method = "L-BFGS-B", 
##     lower = c(0.01, 0.01))
## 
## Coefficients:
##        Estimate  Std. Error
## theta 0.5319061 0.062696418
## sigma 0.1982440 0.004432504
## 
## -2 log L: 120.4586

#LogLik Ornstein-Uhlenbeck
OU.lik <- function(theta1, theta2, theta3){
  n <- length(X)
  dt <- deltat(X)
  - sum(dcOU(X[2:n], dt , X [1:(n -1)],
             c(theta1, theta2, theta3), log = TRUE))
}

set.seed (123)
X <- sde.sim(model = "OU", theta = c(3 ,1 ,2), N = 1000 , delta = 0.01)
autoplot(X, main = 'Simulação para o modelo Ornstein-Uhlenbeck')
fit2 <- mle(OU.lik, start = list(theta1 = 1 , theta2 = 0.5 , theta3 = 1),
            method = "L-BFGS-B", 
            lower = c(-Inf, 0, 0))
summary(fit2)

## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = OU.lik, start = list(theta1 = 1, theta2 = 0.5, 
##     theta3 = 1), method = "L-BFGS-B", lower = c(-Inf, 0, 0))
## 
## Coefficients:
##        Estimate Std. Error
## theta1 5.758780 1.88976587
## theta2 1.778863 0.56974127
## theta3 1.988226 0.04481315
## 
## -2 log L: -410.5437

#LogLik Cox-Ingersoll-Ross
CIR.lik<-function ( theta1 , theta2 , theta3 ) {
  n <-length(X)
  dt <- deltat(X)
  -sum(sde::dcCIR(x = X[2:n], Dt = dt, x0=X[1:(n-1)], 
                  theta =c(theta1, theta2, theta3), log = TRUE))
}

set.seed(123)
X<-sde.sim(X0 =.1, model="CIR", 
             theta = c(.2 , 0.06 , 0.15), N = 1000 , delta = 0.01)
autoplot(X, main = 'Simulação para o modelo Cox-Ingersoll-Ross')
fit3<-mle(CIR.lik , start = list(theta1 =.1 , theta2 =.1, theta3 =.3),
      method = "L-BFGS-B" , lower = c(0.001 ,0.001 ,0.001),
      upper =c (1 ,1 ,1))
summary(fit3)

## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = CIR.lik, start = list(theta1 = 0.1, theta2 = 0.1, 
##     theta3 = 0.3), method = "L-BFGS-B", lower = c(0.001, 0.001, 
##     0.001), upper = c(1, 1, 1))
## 
## Coefficients:
##         Estimate  Std. Error
## theta1 0.1705107 0.056418298
## theta2 0.0861770 0.098528222
## theta3 0.1423016 0.003182266
## 
## -2 log L: -6207.915

## ❯ Scraping historical crypto data
## ❯ Processing historical crypto data