1 VaR por simulación de Monte carlo.
- 1.1 Proceso estocástico para el precio de una acción.
2 Distribución del ultimo día simulado.
3 Convergencia de la distribución.
4 Convergencia de la distribución.
5 Aplicación.
6 generación trazas de precios correlacionadas.
7 Modelos dinámicos para la estimación de VaR.

1 VaR por simulación de Monte carlo.

La estimación del VaR por simulación de Monte carlo es una metodología que consiste en generar escenarios aleatorios a través de muestreos aleatorios de distribuciones de probabilidad conocidas que se supone que describen el comportamiento de las distribuciones de perdida de los factores de riesgo.
Se requiere generar una cantidad suficiente de escenarios con el fin de que la convergencia de la distribución simulada hacia la distribución real, sea lo suficientemente precisa para la aplicación requerida.
Se puede utilizar cualquier distribución de probabilidad para los factores de riesgo.
Posibilita modelar cualquier portafolio compuesto por instrumentos lineales y no lineales.

1.1 Proceso estocástico para el precio de una acción.

Para la simulación del precio de acciones es común suponer que estas siguen un paseo aleatorio, el cual puede ser representado con un proceso estocástico Browniano geométrico.

\[dS=\mu Sdt+\sigma S dw\]

\(S\) corresponde al precio de una acción o valor de un portafolio de acciones. \(\mu:\) Corresponde con el valor esperado de los retornos y \(\sigma\) es la volatilidad de los retornos.

\(dw\) Es el diferencial de un proceso de Wiener que tiene las siguientes propiedades:

El diferencial \(dw=\epsilon\sqrt{dt}\), donde \(\epsilon \sim \mathcal{N}(0,1)\).
Los cambios \(dw\) son independientes para cualquier par de intervalos diferentes del tiempo¹

En términos discretos, el delta del proceso browniano geométrico (o \(it\hat{o}\)) se puede escribir así:

\[\Delta S=\mu S \Delta t +\sigma S\epsilon\sqrt{\Delta t}\]

\[S_{t+1}=S_t* e^{ \left[ \left(\mu-\frac{\sigma^2}{2} \right)\Delta t+\sigma\epsilon_t\sqrt{\Delta t} \right]}\]

Ejemplo:

So=100; volatilidad=0.15 anual; mu=0.50 anual ;delta_t=1/365 ;iteraciones=1000; periodos=365

So=100; volatilidad=0.15;  mu=0.50 ;delta_t=1/365 ;iteraciones=1000; periodos=365

matriz_precios=matrix(,periodos+1,iteraciones) # en esta matriz se almacenan las trazas de precios que se generan en la simulación. Cada fila es un momento del tiempo y cada columna es una traza simulada. 

matriz_precios[1,]=So # en la primera fila de todas las columnas, se pone el precio inicial de la simulación.

for(i in 1:iteraciones)  # este for va iterando de columna en columna
for (j in 2:nrow(matriz_precios)) # este for itera entre cada fila. La fila 1 ya contiene el precio inicial.
  {

# Como la volatilidad y el valor esperado de los rendimientos están expresados de manera anual y queremos caminar en pasos de un dia, al multiplicar mu por 1/365 se convierte la tasa anual a tasa diaria y al multiplicar la volatildad por raiz(1/365) la volatilidad anual se escala a volatilidd diaria. 
  

matriz_precios[j,i]=matriz_precios[j-1,i]*exp((mu-volatilidad^2/2)*delta_t+volatilidad*rnorm(1)*sqrt(delta_t))
    
}

drift=apply(matriz_precios, 1, mean)

cuantiles=t(apply(matriz_precios, 1,quantile, probs=c(0.01,0.05,0.95,0.99))) 

matplot(matriz_precios,type="l", col="gray",ylab="")
lines(drift,type="l",lwd=2)
lines(cuantiles[,"1%"],col="blue")
lines(cuantiles[,"99%"],col="blue")
lines(cuantiles[,"5%"],col="red")
lines(cuantiles[,"95%"],col="red")
legend("topleft",legend = c("drift","intervalo 90%","intervalo 98%"),lwd = c(2,1,1,1),bty="n",col=c("black","red","blue"))

library(knitr)

2 Distribución del ultimo día simulado.

Cada una de las filas de la matriz de precios simuladas, es una muestra simulada de los precios para el final de un periodo determinado. Como ejemplo, tomaremos la ultima fila que corresponde a los precios simulados un año adelante (luego de haber caminado 365 días).

dist_Rtos_simulada=log(matriz_precios[nrow(matriz_precios),]/matriz_precios[1,]) # se estiman los rendimientos para el último día como el logaritmo natural del precio simulado para dentro de un año, con el precio incial de la simulación (última observación de la base de datos histórica).


par(mfrow=c(1,2))

hist(dist_Rtos_simulada,100,col="black",main="rtos día 365",xlab="") 
hist(matriz_precios[nrow(matriz_precios),],breaks = 100,col="black",main="precio día 365", xlab = "")

#mu-volatilidad^2/2

Se supone que la distribución de precios simulados corresponde a una distribución log normal y la de los rendimientos a una distribución normal. El VaR a una año se puede calcular como un cuantil empírico de la distribución de retornos simulada.

3 Convergencia de la distribución.

Entre mayor número de iteraciones, la distribución futura de los precios y los rendimientos convergen a su distribución real. por ejemplo si se supone que la distribución futura real de los rendimientos es normal, entonces se puede expresar así:

\[ln(S_T)-ln(S_0)\sim \mathcal{N}(\left(\mu-\frac{\sigma^2}{2} \right)\Delta t,\sigma^2T)\]

4 Convergencia de la distribución.

Así, dadas las condiciones del ejemplo, la distribución de probabilidad para el dia 365 tiene los siguientes parámetros:

\[E(\Delta ln(S_T))=0.5-0.15^2/2=0.48875\]

\[\ \sigma_{r_{T}}=0.15*\sqrt{365/365}=0.15\]

y el cuantil 1% suponiendo normalidad (\(z_{1\%}=-2.326\)) es: \(-2.326*0.15+0.48875=0.1398\)

Tabla de parámetros obtenidos en la simulación

# print("media")
# print(mean(dist_Rtos_simulada))
# 
# print("volatilidad")
# print(sd(dist_Rtos_simulada))
# 
# print("quantiles")

#quantile(dist_Rtos_simulada,0.01)

tabla_param_dist_Rtos_simulada=rbind(c("media","desv","1%"),round(c(mean(dist_Rtos_simulada),sd(dist_Rtos_simulada),quantile(dist_Rtos_simulada,c(0.01))),4))

colnames(tabla_param_dist_Rtos_simulada)=c("","","")

library(knitr)

(tabla_param_dist_Rtos_simulada)

##                                
## [1,] "media"  "desv"   "1%"    
## [2,] "0.4878" "0.1478" "0.1525"

A medida que se aumenta el número de iteraciones, el parámetro converge con mayor exactitud al valor teórico

So=100; volatilidad=0.15;  mu=0.50 ;delta_t=1 ;iteraciones=1000; periodos=1



percentil_1=vector()

#vector_iteraciones=c(seq(10:1000),5000,10000)

for (k in seq(10:1000) ) {

  matriz_precios_2=matrix(,periodos+1,k)
matriz_precios_2[1,]=So  

for(i in 1:k){
for (j in 2:nrow(matriz_precios_2)) {
  
matriz_precios_2[j,i]=matriz_precios_2[j-1,i]*exp((mu-volatilidad^2/2)*delta_t+volatilidad*rnorm(1)*sqrt(delta_t))
    
}}

rtos=diff(log(matriz_precios_2))

percentil_1[k]=quantile(rtos,0.01)

}

percentil_1=percentil_1[!is.na(percentil_1)]

plot(percentil_1,type="l", xlab = "iteraciones",main="percentil 1  en función del número de iteraciones")

El percentil 1 converge a su valor teórico con mayor precisión cuando se aumenta el número de iteraciones.

5 Aplicación.

En este documento se hace la simulación de Montecarlo de la distribución de perdidas de un portafolio de inversión en acciones. Se utiliza una base de precios diarios de un año desde octubre 1 de 2019 hasta septiembre 30 de 2020.

Las acciones analizadas son: “ECOPETROL”, “PFBCOLOM” , “ISA”, “GRUPOSURA” ,“GEB” “NUTRESA”.

La simulación de Monte carlo multivariada, requiere que la estructura de correlaciones entre las variables aleatorias generadas, tengan la estructura de correlaciones que se supone que hay entre los factores de riesgo del portafolio.

base_datos_portafolio <- read.csv("~/Library/Mobile Documents/com~apple~CloudDocs/curso_riesgos_-msc_finanzas_R/base_datos_portafolio.csv", sep=";")


library(matrixStats)

Rtos_acciones=colDiffs(as.matrix(log(base_datos_portafolio[,-c(1,2)])))  


colnames(Rtos_acciones)=colnames(base_datos_portafolio)[3:8]

5.1 Gráficas de precios y rendimientos.

La series presentan un salto en marzo de 2020 debido al impacto del COVID-19.

par(mfrow=c(3,2))


for(j in 3:ncol(base_datos_portafolio)){
  
plot(base_datos_portafolio[,j],type="l",main=names(base_datos_portafolio)[j],ylab="")
  
}

par(mfrow=c(3,2))

for (i in 1:ncol(Rtos_acciones)) {
  

plot(Rtos_acciones[,i],type="l",ylab="", main = colnames(Rtos_acciones)[i])


}

5.2 matriz de correlaciones.

matriz_correlaciones=cor(Rtos_acciones)



kable(round(matriz_correlaciones,3))

	ECOPETROL	PFBCOLOM	ISA	GRUPOSURA	GEB	NUTRESA
ECOPETROL	1.000	0.600	0.453	0.622	0.384	0.517
PFBCOLOM	0.600	1.000	0.573	0.734	0.496	0.497
ISA	0.453	0.573	1.000	0.649	0.523	0.479
GRUPOSURA	0.622	0.734	0.649	1.000	0.489	0.567
GEB	0.384	0.496	0.523	0.489	1.000	0.398
NUTRESA	0.517	0.497	0.479	0.567	0.398	1.000

5.3 Descomposición de Cholesky.

La factorización de Cholesky de una matriz de correlaciones \(\rho\), consiste en expresarla como producto de una matriz L triangular por su transpuesta.

\[\varrho=LL^T\]

En el caso de 2 activos:

\[L=\begin{bmatrix}1 & \rho_{12} \\ 0 & \sqrt{1-\rho_{12}^{2}} \end{bmatrix}\]

5.4 Matriz triangular L

L=chol(matriz_correlaciones)

kable(round(L,3))

	ECOPETROL	PFBCOLOM	ISA	GRUPOSURA	GEB	NUTRESA
ECOPETROL	1	0.6	0.453	0.622	0.384	0.517
PFBCOLOM	0	0.8	0.376	0.451	0.332	0.234
ISA	0	0.0	0.808	0.244	0.276	0.194
GRUPOSURA	0	0.0	0.000	0.591	0.054	0.156
GEB	0	0.0	0.000	0.000	0.814	0.073
NUTRESA	0	0.0	0.000	0.000	0.000	0.782

Si \(Z_1^t\) y \(Z_2^t\) son dos variables aleatorias independientes generadas para el periodo \(t\) entonces:

\[(Z_1^t,Z_2^t)\begin{bmatrix}1 & \rho_{12} \\ 0 & \sqrt{1-\rho_{12}^{2}} \end{bmatrix}=(Z_1^t,\rho_{12}Z_1^t+\sqrt{1-\rho_{12}^{2}}Z_2^t)=(k_1,k_2)\]

6 generación trazas de precios correlacionadas.

Para el ejemplo se utilizan los siguientes pesos:

ECOPETROL	PFBCOLOM	ISA	GRUPOSURA	GEB	NUTRESA
0.148	0.184	0.120	0.113	0.156	0.279

Las volatilidades, el valor esperado y las correlaciones se estiman de manera histórica, utilizando todo el histórico disponible en la base de datos.

W=c(0.148, 0.184, 0.120, 0.113, 0.156, 0.279)

Volatilidades=apply(Rtos_acciones,2, sd) # vector que contiene las volatilidades históricas de cada acción.

precios_iniciales=base_datos_portafolio[nrow(base_datos_portafolio),-c(1,2)] #Vector que contiene los precios iniciales de cada acción 

media_rtos=apply(Rtos_acciones,2, mean) # vector que contiene la media histórica de cada acción.

delta_t=10 #  La volatilidad es diaria y queremos estimar el VaR para dentro de 10 dias y solo vamos a caminar un paso de 10 días.

iteraciones=1000; 


precio_activos=matrix(,2,length(Volatilidades)) # en esta matriz se almacenan temporalmente los precios iniciales (fila 1) y los precios finales que corresponden con el día 10 (fila 2).

 
precio_activos[1,]=t(as.vector(precios_iniciales)) # se asigna el precio inicial de cada acción.

rtos_portafolio=vector() # En este vector se almacenan los rendimientos simulados del portafolio. 

j=1
i=1
for(j in 1:iteraciones)  {
for (i in 1:length(Volatilidades))   {
  
Z=rnorm(length(Volatilidades)) # vector de variables aleatorios estándar independientes

K=Z%*%L  # vector de variables aleatorios  correlacionadas.

precio_activos[2,i]=precio_activos[1,i]*exp((media_rtos[i]-Volatilidades[i]^2/2)*delta_t+Volatilidades[i]*K[i]*sqrt(delta_t)) # ecuación de paseo aleatorio. 


Rtos_activos=colDiffs(log(precio_activos)) 

}
rtos_portafolio[j]=(Rtos_activos*W)  
  }  


VaR_10_99=quantile(-rtos_portafolio,0.99) # La función quantile estima el cuantil empirico de una distribución de frecuencias.


print("vaR 10 días  del portafolio 99% de confianza")

## [1] "vaR 10 días  del portafolio 99% de confianza"

print(VaR_10_99)

##        99% 
## 0.03781581

hist(-rtos_portafolio,100,col="black", main="distribución de perdidas y ganancias del portafolio" )
abline(v=VaR_10_99, col="blue")
text(VaR_10_99, 20,round(VaR_10_99,4) ,cex = 0.65)

EL VaR se estima como el 1-\(\alpha\) cuantil de la distribución de frecuencias simulada .

7 Modelos dinámicos para la estimación de VaR.

Los modelos de media y volatilidad dinámica, permiten hacer estimaciones mas ajustadas a las condiciones actuales de una serie de tiempo. En el caso de que existan autocorrelaciones entre los rendimientos y autocorrelaciones del cuadrado de los rendimientos, es posible estimar estos parámetros en función de rezagos.

7.1 Modelo RiskMetrics.

Modelo de varianza:

Longerstaey and Spencer (1996)

\[\sigma_{t+1}^{2}=\lambda \sigma_{t}^{2}+(1-\lambda)r_{t}^{2}\]

Bajo la metodología de Riskmetrics se estima el VaR bajo el supuesto de normalidad, suponiendo que \(\mu=0\) y \(\sigma\) constante, se reemplaza por \(\sigma_t\) dinámico.

Para un horizonte de \(T\) días \(\sigma_{t+T}^{2}=T\sigma_{t+1}^{2}\)

\[VaR_{h,\alpha}\approx \Phi^{-1}(1-\alpha)\sqrt{T}\sigma_t\]

\[\sigma_{t}^{2}=\lambda \sigma_{t-1}^{2}+(1-\lambda)r_{t-1}^{2}\]

7.2 Aplicación:

Se aplica el modelo sobre los rendimientos de las acciones del portafolio que se viene tratando en este documento. Se utliza arbitrariamente un \(\lambda\) de 0.9 y se sigue suponiendo que las correlaciones permanecen constantes. El VaR se estima con el mismo procedimiento aplicado con a estimación histórica de la volatilidad, pero en este caso se usa la volatilidad dinámica \(\sigma_{t}\)

lambda_var=0.9

Rtos_acciones_cuadrado=Rtos_acciones^2

varianza_risk_metric=matrix(,nrow(Rtos_acciones_cuadrado),ncol(Rtos_acciones_cuadrado)) # en esta matriz se almacena la volatilidad estimada.

varianza_risk_metric[1,]=Volatilidades^2*lambda_var+(1-lambda_var)*Rtos_acciones_cuadrado[1,] # se inicia la serie de estimación de la varianza con la varianza histórica de la muestra.

### se estiman las volatilidades para el resto de periodos.

for (j in 2:nrow(Rtos_acciones_cuadrado)) {

varianza_risk_metric[j,]=varianza_risk_metric[j-1,]*lambda_var+(1-lambda_var)*Rtos_acciones_cuadrado[j,] 

sigma_risk_metric=sqrt(varianza_risk_metric)


}

## VaR riskmetrics.


Volatilidades_risk_m=sigma_risk_metric[nrow(sigma_risk_metric),] # se selecciona la volatilidad estimada para un periodo adelante del ultimo dato observado para hacer las proyecciones con la simulación montecarlo

precio_activos_riskM=matrix(,2,length(Volatilidades_risk_m)) # en esta matriz se almacenan temporalmente los precios iniciales (fila 1) y los precios finales que corresponden con el día 10 (fila 2).

 
precio_activos_riskM[1,]=t(as.vector(precios_iniciales)) # se asigna el precio inicial de cada acción a la fila 1 de la matriz de precios.

rtos_portafolio_risk_m=vector() # En este vector se almacenan los rendimientos simulados del portafolio. 

for(j in 1:iteraciones)  {
for (i in 1:length(Volatilidades_risk_m))   {
  
Z=rnorm(length(Volatilidades_risk_m)) # vector de variables aleatorios estándar independientes

K=Z%*%L  # vector de variables aleatorios  correlacionadas.

precio_activos_riskM[2,i]=precio_activos_riskM[1,i]*exp((media_rtos[i]-Volatilidades_risk_m[i]^2/2)*delta_t+Volatilidades_risk_m[i]*K[i]*sqrt(delta_t)) # ecuación de paseo aleatorio. 


Rtos_activos=colDiffs(log(precio_activos_riskM)) # vetor de rendimientos calculado como el el logaritmo del precio simuado para el dia 10 y el precio inicial de la simulación (que corresponde con el último observado en la muestra)

}
rtos_portafolio_risk_m[j]=(Rtos_activos*W) # se  calcula la rentabilidad simulada del portafolio para la traza j. El vector rtos_portafolio_risk_m contiene cada uno de los rendimientos simulados en cada iteración.   
  }  


VaR_10_99_rm=quantile(-rtos_portafolio_risk_m,0.99) # La función quantile estima el cuantil empirico de una distribución de frecuencias.


print("vaR risk metrics 10 días  del portafolio 99% de confianza")

## [1] "vaR risk metrics 10 días  del portafolio 99% de confianza"

print(VaR_10_99_rm)

##        99% 
## 0.02594589

hist(-rtos_portafolio_risk_m,100,col="black", main="distribución de perdidas y ganancias del portafolio" )
abline(v=VaR_10_99_rm, col="blue")
text(VaR_10_99_rm, 20,round(VaR_10_99_rm,4) ,cex = 0.65)

Volatilidades de cada activo estimadas con Risk Metrics

matplot(sigma_risk_metric, type="l", main="Volatilidades Risk Metrics")

7.3 Modelo Riskmetrics multivariado.

La metodología RiskMetrics construye pronosticos de covarianza y correlación de la misma manera que lo hace con el pronostico de volatilidad volatilidad, pero en vez de trabajar con una serie de retornos al cuadrado, trabaja con el producto de 2 series de retornos, de tal manera que la covarianza entre un factor \(i\) y un factor \(j\) estimada en \(t\) para un periódo \(t+1\) es:

\[\sigma_{i,j,t+1}=\lambda \sigma_{i,j,t}+(1-\lambda)r_{i,t}r_{j,t}\]

Escalamiento de la covarianza para un horizonte de \(T\) días.

\[\sigma_{i,j,t+T}=T\sigma_{i,j,t+1}\]

7.4 Correlación Riskmetrics

\[\rho_{i,j,t+1}=\frac{\sigma_{i,j,t+1}}{\sigma_{i,t+1}.\sigma_{j,t+1}}\]

para un horizonte T:

\[\rho_{i,j,t+T}=\frac{T\sigma_{i,j,t+1}}{\sqrt{T}\sigma_{i,t+1}.\sqrt{T}\sigma_{j,t+1}}=\rho_{i,j,t+1}\]

La matriz de correlaciones permanece invariante en el tiempo, independiente del horizonte de estimación.

En terminos matriciales la matriz de correlaciones en el momento \(t\) se estima así:

\[\varrho_t=S_t^{-1}\Sigma_t S_t^{-1}\]

\[\Sigma_t=\begin{bmatrix}\sigma_{t1}^2 & \sigma_{t1,2} &. & \sigma_{t1,j} &. & \sigma_{t1,n}\\\sigma_{t2,1} & \sigma_{t2}^2 &\sigma_{t2,j} &. &.& \sigma_{t2,n} \\. &. &. &. &. &.\\\sigma_{tj,1} & \sigma_{tj,2} &. & \sigma_{tj}^2 &. &\sigma_{tj,n}\\. &. &. &. &. &. \\\sigma_{tn,1} & \sigma_{tn,2} & \sigma_{tn,j} &. &. & \sigma_{tn}^2 \end{bmatrix} \]

Es la matriz de covarianzas.

\(S_t\) es una matriz diagonal:

\[S_t=\begin{bmatrix}\sigma_{t1} & 0 &. & 0 &. & 0\\0 & \sigma_{t2} &. &0 &.& 0\\. &. &. &. &. &. \\0 & 0 & &. \sigma_{tj} &. &0\\. &. &. &. &. &. \\0 & 0 & 0 &. & &. \sigma_{tn} \end{bmatrix}; \]

\[\varrho_t=\begin{bmatrix}1 & \rho_{t1,2} &. & \rho_{t1,j} &. & \rho_{t1,n}\\\rho_{t2,1} & 1 &\rho_{t2,j} &. &.& \rho_{t2,n} \\. &. &. &. &. &.\\\rho_{tj,1} & \rho_{tj,2} &. & 1 &. &\rho_{tj,n}\\. &. &. &. &. &. \\\rho_{tn,1} & \rho_{tn,2} & \rho_{tn,j} &. &. & 1 \end{bmatrix} \]

library(RiskPortfolios) #esta libreria estima la varianza yla covarianza EWMA


matriz_cov_risk_m=covEstimation(Rtos_acciones, control = list(type = 'ewma', lambda = 0.9))

# covEstimation es una función de la libreria RiskPortfolios que estima la matriz de covarianzas  EWMA

S=sqrt(diag(diag(matriz_cov_risk_m))) #S es una matriz diagonal que contiene las desviaciones estandar de cada activo del portafolio.


matriz_correlaciones_risk_m= solve(S)%*%matriz_cov_risk_m%*%solve(S)


L_rism_m=chol(matriz_correlaciones_risk_m)

colnames(matriz_correlaciones_risk_m)=colnames(matriz_correlaciones)

rownames(matriz_correlaciones_risk_m)=rownames(matriz_correlaciones)

kable(matriz_correlaciones_risk_m)

	ECOPETROL	PFBCOLOM	ISA	GRUPOSURA	GEB	NUTRESA
ECOPETROL	1.0000000	0.3944260	-0.0156800	0.3120838	0.3654095	-0.1988378
PFBCOLOM	0.3944260	1.0000000	0.5353386	0.5516625	0.1549129	0.1171439
ISA	-0.0156800	0.5353386	1.0000000	0.4569694	-0.0328785	0.2731609
GRUPOSURA	0.3120838	0.5516625	0.4569694	1.0000000	0.0972055	0.1238754
GEB	0.3654095	0.1549129	-0.0328785	0.0972055	1.0000000	0.1803432
NUTRESA	-0.1988378	0.1171439	0.2731609	0.1238754	0.1803432	1.0000000

Matriz de correlaciones Riskmetrics

kable(matriz_correlaciones)

	ECOPETROL	PFBCOLOM	ISA	GRUPOSURA	GEB	NUTRESA
ECOPETROL	1.0000000	0.5995216	0.4533668	0.6222821	0.3841028	0.5166334
PFBCOLOM	0.5995216	1.0000000	0.5730359	0.7344088	0.4962210	0.4967128
ISA	0.4533668	0.5730359	1.0000000	0.6494365	0.5225948	0.4790865
GRUPOSURA	0.6222821	0.7344088	0.6494365	1.0000000	0.4887676	0.5666444
GEB	0.3841028	0.4962210	0.5225948	0.4887676	1.0000000	0.3975256
NUTRESA	0.5166334	0.4967128	0.4790865	0.5666444	0.3975256	1.0000000

Matriz de correlaciones historicas.

#####

Volatilidades_risk_m_cov=sqrt(diag(matriz_cov_risk_m))

precio_activos_riskM_cov=matrix(,2,length(Volatilidades_risk_m_cov)) # en esta matriz se almacenan temporalmente los precios iniciales (fila 1) y los precios finales que corresponden con el día 10 (fila 2).

 
precio_activos_riskM_cov[1,]=t(as.vector(precios_iniciales)) # se asigna el precio inicial de cada acción.

rtos_portafolio_risk_m_cov=vector() # En este vector se almacenan los rendimientos simulados del portafolio. 

for(j in 1:iteraciones)  {
for (i in 1:length(Volatilidades_risk_m_cov))   {
  
Z=rnorm(length(Volatilidades_risk_m_cov)) # vector de variables aleatorios estándar independientes

K=Z%*%L_rism_m  # vector de variables aleatorios  correlacionadas.

precio_activos_riskM_cov[2,i]=precio_activos_riskM_cov[1,i]*exp((media_rtos[i]-Volatilidades_risk_m_cov[i]^2/2)*delta_t+Volatilidades_risk_m_cov[i]*K[i]*sqrt(delta_t)) # ecuación de paseo aleatorio. 


Rtos_activos=colDiffs(log(precio_activos_riskM_cov)) 

}
rtos_portafolio_risk_m_cov[j]=(Rtos_activos*W)  
  }  


VaR_10_99_rm_cov=quantile(-rtos_portafolio_risk_m_cov,0.99) # La función quantile estima el cuantil empirico de una distribución de frecuencias.


print("vaR risk metrics cov 10 días  del portafolio 99% de confianza")

## [1] "vaR risk metrics cov 10 días  del portafolio 99% de confianza"

print(VaR_10_99_rm_cov)

##        99% 
## 0.02487223

hist(-rtos_portafolio_risk_m_cov,100,col="black", main="distribución de perdidas y ganancias del portafolio" )
abline(v=VaR_10_99_rm_cov, col="blue")
text(VaR_10_99_rm_cov, 20,round(VaR_10_99_rm_cov,4) ,cex = 0.65)

Se observan diferencias importantes en la matriz de correlaciones histórica y la Riscmetrics debido al choque ocasionado por el COVID-19.

7.5 VaR ARMA-Garch

Dada la naturaleza autorregresiva de la media y la varianza de losfactores de riesgo financieros, es posible estimar la media y lavarianza de las distribuciones de perdidas a través de modelosARMA y GARCH respectivamente.

library(stats)
library(zoo)

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(graphics)
library(quadprog)
library(tseries)

## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

acf(Rtos_acciones[,1],main=colnames(Rtos_acciones)[1])

pacf(Rtos_acciones[,1],main=colnames(Rtos_acciones)[1])

Modelo de media

modelomedia_1 <- arma(Rtos_acciones[,1], lag=list(ar=c(1)),include.intercept = FALSE)

## Warning in arma(Rtos_acciones[, 1], lag = list(ar = c(1)),
## include.intercept = FALSE): order is ignored

## Warning in optim(coef, err, gr = NULL, hessian = TRUE, ...): one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly

summary(modelomedia_1)

## 
## Call:
## arma(x = Rtos_acciones[, 1], lag = list(ar = c(1)), include.intercept = FALSE)
## 
## Model:
## ARMA(1,0)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.2075159 -0.0111579  0.0007874  0.0123848  0.1303260 
## 
## Coefficient(s):
##      Estimate  Std. Error  t value Pr(>|t|)    
## ar1   0.24070     0.06225    3.866  0.00011 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Fit:
## sigma^2 estimated as 0.001122,  Conditional Sum-of-Squares = 0.27,  AIC = -959.03

modelomedia_1_residual=modelomedia_1$residuals

modelomedia_1_residual=modelomedia_1_residual[!is.na(modelomedia_1_residual)]

acf(modelomedia_1_residual)

acf(modelomedia_1_residual^2)

Modelo de volatilidad

modelovarianza_1<-garch(modelomedia_1_residual, order=c(0,1) )

## 
##  ***** ESTIMATION WITH ANALYTICAL GRADIENT ***** 
## 
## 
##      I     INITIAL X(I)        D(I)
## 
##      1     1.065787e-03     1.000e+00
##      2     5.000000e-02     1.000e+00
## 
##     IT   NF      F         RELDF    PRELDF    RELDX   STPPAR   D*STEP   NPRELDF
##      0    1 -7.063e+02
##      1    6 -7.072e+02  1.25e-03  1.79e-03  1.0e-03  1.3e+08  1.0e-04  1.13e+05
##      2    7 -7.073e+02  1.48e-04  2.04e-04  8.2e-04  2.0e+00  1.0e-04  2.01e+01
##      3    8 -7.073e+02  3.87e-05  4.04e-05  9.9e-04  2.0e+00  1.0e-04  1.98e+01
##      4   15 -7.210e+02  1.90e-02  2.38e-02  5.1e-01  2.0e+00  1.0e-01  1.94e+01
##      5   16 -7.304e+02  1.28e-02  1.25e-02  2.5e-01  6.5e-01  1.0e-01  1.59e-02
##      6   17 -7.387e+02  1.12e-02  1.28e-02  1.7e-01  7.8e-01  1.0e-01  1.57e-02
##      7   19 -7.449e+02  8.41e-03  6.51e-03  2.0e-01  0.0e+00  1.8e-01  6.51e-03
##      8   20 -7.474e+02  3.27e-03  2.21e-03  9.8e-02  0.0e+00  1.2e-01  2.21e-03
##      9   21 -7.485e+02  1.58e-03  1.25e-03  8.2e-02  1.6e-01  1.2e-01  1.26e-03
##     10   22 -7.489e+02  5.31e-04  4.14e-04  5.8e-02  0.0e+00  9.4e-02  4.14e-04
##     11   23 -7.490e+02  7.56e-05  6.35e-05  2.6e-02  0.0e+00  4.6e-02  6.35e-05
##     12   24 -7.490e+02  4.05e-06  3.86e-06  6.5e-03  0.0e+00  1.2e-02  3.86e-06
##     13   25 -7.490e+02  5.81e-08  7.92e-08  8.4e-04  0.0e+00  1.5e-03  7.92e-08
##     14   26 -7.490e+02  5.18e-09  5.20e-09  1.0e-04  0.0e+00  1.9e-04  5.20e-09
##     15   27 -7.490e+02  4.54e-13  4.59e-13  1.6e-06  0.0e+00  2.9e-06  4.59e-13
## 
##  ***** RELATIVE FUNCTION CONVERGENCE *****
## 
##  FUNCTION    -7.490054e+02   RELDX        1.592e-06
##  FUNC. EVALS      27         GRAD. EVALS      16
##  PRELDF       4.595e-13      NPRELDF      4.595e-13
## 
##      I      FINAL X(I)        D(I)          G(I)
## 
##      1    3.693359e-04     1.000e+00    -2.787e-03
##      2    9.223979e-01     1.000e+00    -2.170e-06

parametros_modelovarianza=summary(modelovarianza_1)

parametros_modelovarianza$coef

##        Estimate   Std. Error   t value     Pr(>|t|)
## a0 0.0003693359 3.249832e-05 11.364770 0.000000e+00
## a1 0.9223979474 1.246912e-01  7.397456 1.387779e-13

media_estimada=modelomedia_1$fitted.values[-c(1,2)]
volatilidad_estimada=modelovarianza_1$fitted.values[,1]

volatilidad_estimada=volatilidad_estimada[!is.na(volatilidad_estimada)]

VaR_1_99=qnorm(0.01)*volatilidad_estimada+media_estimada


plot(Rtos_acciones[-c(1,2),1],type="l",ylim=c(-0.50,0.11))
lines(VaR_1_99,type="l",col="blue")
legend("topleft",legend = c("rentabilidad observada","VaR_1dia_99%"),col=c("black","blue"),lty=c(1,1),bty = "n")

Longerstaey, Jacques, and Martin Spencer. 1996. “Riskmetricstm—Technical Document.” Morgan Guaranty Trust Company of New York: New York 51: 54.

la independencia lineal implica que: \(cov(dw_i,dw_j=0)\), aunque el supuesto implica independencia en cualquier sentido↩

simulación Monte carlo portafolio

juan Fernando Rendón

10/19/2020