Librerías

library(moments)
library(nortest)
library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(quantmod)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## Loading required package: TTR
library(zoo)
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ tibble  3.1.6     ✓ purrr   0.3.4
## ✓ tidyr   1.2.0     ✓ stringr 1.4.0
## ✓ readr   2.1.1     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x lubridate::as.difftime() masks base::as.difftime()
## x lubridate::date()        masks base::date()
## x dplyr::filter()          masks stats::filter()
## x xts::first()             masks dplyr::first()
## x lubridate::intersect()   masks base::intersect()
## x dplyr::lag()             masks stats::lag()
## x xts::last()              masks dplyr::last()
## x lubridate::setdiff()     masks base::setdiff()
## x lubridate::union()       masks base::union()
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:purrr':
## 
##     some
## The following object is masked from 'package:dplyr':
## 
##     recode
library(astsa)
library(foreign)
library(lmtest)
library(dynlm)
library(nlme)
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
library(lmtest)
library(broom)
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(knitr)
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
library(parallel)
library(car)
library(mlogit)
## Loading required package: dfidx
## 
## Attaching package: 'dfidx'
## The following object is masked from 'package:MASS':
## 
##     select
## The following object is masked from 'package:stats':
## 
##     filter
library(dplyr)
library(tidyr)
library(forecast)
## 
## Attaching package: 'forecast'
## The following object is masked from 'package:nlme':
## 
##     getResponse
## The following object is masked from 'package:astsa':
## 
##     gas
library(fpp2)
## ── Attaching packages ────────────────────────────────────────────── fpp2 2.4 ──
## ✓ fma       2.4     ✓ expsmooth 2.3
## ── Conflicts ───────────────────────────────────────────────── fpp2_conflicts ──
## x forecast::getResponse() masks nlme::getResponse()
## x car::some()             masks purrr::some()
## 
## Attaching package: 'fpp2'
## The following object is masked from 'package:astsa':
## 
##     oil
library(stats)

Simulación Monte Carlo y Bootstrap

Objetivo:

Aplicar el método Monte Carlo y Bootstrap para realizar el pronóstico de una serie de datos empíricos.

Actividad:

En esta actividad se utilizará la serie seleccionada en la Tarea 2, para aplicar el método

Monte Carlo (Distribución Normal y NIG) así como método Bootstrap para realizar para realizar el pronóstico de una serie de datos empíricos para un periodo de 20 observaciones realizando 10000 realizaciones (simulaciones).

#Obtener información en Yahoo

Televisa<-getSymbols("TV", from="2020-01-01", src = "yahoo", 
                    auto.assign = F)[,6]
str(Televisa)
## An 'xts' object on 2020-01-02/2022-07-01 containing:
##   Data: num [1:630, 1] 11.8 11.7 11.6 11.3 11.5 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr "TV.Adjusted"
##   Indexed by objects of class: [Date] TZ: UTC
##   xts Attributes:  
## List of 2
##  $ src    : chr "yahoo"
##  $ updated: POSIXct[1:1], format: "2022-07-03 03:42:38"

1. Parámetros obtenidos al ajustar la distribución que se utilizarán para las simulaciones, no aplica para el caso del método bootstrap

#se divide la serie en dos subconjuntos, los primeros 610 valores se utilizarán para realizar el forecast y los últimos 20 servirán para verificar la precisión de las simulaciones.
TV1 <- head(Televisa, n = 6100)
#Se obtienen los últimos 50 valores del la serie
TV.observed <- tail(Televisa, n = 20)
#Se obtienen las diferencias lograrítmicas de la series de ajuste.
l1.tv <- diff(log(TV1))
l1.tv.df <- as.data.frame(l1.tv)
#se estima la prueba ADF para rechazar H0: se tiene raíz unitaria.
l1.tv<-na.omit(l1.tv)
tv.adf <- adf.test(l1.tv)
## Warning in adf.test(l1.tv): p-value smaller than printed p-value
tv.adf
## 
##  Augmented Dickey-Fuller Test
## 
## data:  l1.tv
## Dickey-Fuller = -8.151, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary

Como el p-value <-05, entonces rechazamos la H0, es decir, que nuestra serie es estacionaria.

****——————————————————**** #Es importante mencionar que para realizar el ajuste de series de rendimientos debemos considerar distribuciones continuas que su dominio sea R, ya que los rendimeintos se asume tiene media cero, por lo que también se observan n rendimientos negativos. ****——————————————————****

#————-Metodo Monte Carlo (Normal dist.)——————

N = 10000 #número de realizaciones, estados de la naturaleza, trayectorias.
t = 20 #Número de periodos en el futuro
min.Price = 8.75
#forma de aprovechar lamanera en que R crea las matrices
#Se estiman los parámetros de la distribución Normal.
media <- mean(l1.tv)
des <- sd(l1.tv)
#Se crea una matriz en la que se almacenan los valores simulados de la v.a. asumiendo independencia y estacionariedad. Cuando se crea una matriz se declaran tres parámetros 
#1. Valores que contendrá la matriz. 
#2. Número de filas.
#3. Número de columnas.
#En este caso aprovechamos para realizar las simulaciones utilizando la función rnorm() al mismo momento que declaramos la matriz, ya que estableceremos que los valores con los que se debe llenar la matriz son las simulaciones de la v.a.

MC.Norm <- matrix(rnorm(N*t,media,des),t,N)
#se obtiene el valor de los rendimientos acumulados, con lo que se forman las trayectorias de los rendimientos.
MC.Norm.Tra <- apply(MC.Norm,2,cumsum)
#Se obtiene el último precio de la serie (en este caso el valor 610 de la serie)
P0 <- as.numeric(last(TV1))
#Se genera un vector N con el último precio de la serie.
P0.vector <- rep(P0,N)
#debido a que los rendimientos se obtivieron mediante la diferencia logarítmica manera de reescalar al precio es por medio de la exponencial
MC.Norm.Price <- P0*exp(MC.Norm.Tra)
#Se agrega en la primera fila de la matriz un vector con P0 para tener la referencia de inicio de las simulaciones.
MC.Norm.Price <- rbind(P0.vector,MC.Norm.Price)
#Se agrega una columna con el número de periodos de tiempo en los que se observará el precio 
MC.Norm.Price <- cbind(Time=1:(t+1),MC.Norm.Price)
#Se  reornenan los datos para graficarlos
library(tidyverse)
MC.Norm.Price.DF <- as.data.frame(MC.Norm.Price)
MC.Norm.Price.tidy <- MC.Norm.Price.DF %>% 
  gather(-Time,key = Sim,value = Price)

MC.Norm.Price.tidy %>% 
  ggplot(aes(x=Time,
             y=Price,
             group=Sim))+
  geom_line(alpha=0.1,
            color="blue")+
  theme(legend.position = 'none')+
  labs(title = 'MC Normal Simulation',
       subtitle = 'Televisa/USD',
       x="Days",
       y="Price")+
  geom_hline(yintercept = min.Price,
             col="red",
             size=.2)+
  theme_bw()

#Estimación del Valor Esperado
#Se crea el objeto valor esperado
EV.Norm = NA
#Se estima el valor esperado para cada t
for(i in 1:(t+1)){
  EV.Norm[i] <- sum(MC.Norm.Price.DF[i,-1])/N
}
#Se crea un Data Frame con los Valores esperados estimados
EV.Norm.DF <- tibble(Time = 1:(t+1), Esperado =  EV.Norm)
#Se genera nuevamente la gráfica de trayectorias agregando una linea de Valor Esperado
MC.Norm.Price.tidy %>% 
  ggplot(aes(x=Time,
             y=Price,
             group=Sim))+
  geom_line(alpha=0.1,
            color="blue")+
  theme(legend.position = 'none')+
  labs(title = 'MC Normal Simulation',
       subtitle = 'Televisa/USD',
       x="Days",
       y="Price")+
  geom_hline(yintercept = min.Price,
             col="red",
             size=.2)+
  geom_line(data = EV.Norm.DF,inherit.aes = F, aes(x= Time, y=Esperado), size =1)+
  theme_bw()

#Estimación de quantiles para un nivel de confianza del 95%
MC.Norm.Qs = NA
MC.Norm.Qi = NA

for(i in 1:(t+1)){
  MC.Norm.Qi[i] <- as.numeric(quantile(MC.Norm.Price.DF[i,-1], c(0.025,0.975))[1])
  MC.Norm.Qs[i] <- as.numeric(quantile(MC.Norm.Price.DF[i,-1], c(0.025,0.975))[2])
}
## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames
MC.Norm.BandasQ <- tibble(Time =1:(t+1), sup = MC.Norm.Qs, inf = MC.Norm.Qi)
MC.Norm.Price.tidy %>% 
  ggplot(aes(x=Time,
             y=Price,
             group=Sim))+
  geom_line(alpha=0.1,
            color="blue")+
  theme(legend.position = 'none')+
  labs(title = 'MC Normal Simulation',
       subtitle = 'Televisa/USD',
       x="Days",
       y="Price")+
  geom_hline(yintercept = min.Price,
             col="red",
             size=.2)+
  geom_ribbon(data = MC.Norm.BandasQ,inherit.aes = F, aes(x= Time, ymin = inf, ymax = sup), fill = "purple", alpha = 0.6)+
  geom_line(data = EV.Norm.DF,inherit.aes = F, aes(x= Time, y=Esperado), size =1)+
  theme_bw()

#Se crea una matriz donde se almacenarán los resultados del valor esperado para su 
#posterior análisis
Resultados = matrix(0,6,21)
#Se almacenan los valores empíricos en la matriz
Resultados[1,] <- c(P0,TV.observed)
#Se almacenan los Valores Esparados estimados con la Normal
Resultados[2,] <- EV.Norm

————–Metodo Monte Carlo (t-Student)———

N = 10000 #número de realizaciones, estados de la naturaleza, trayectorias.
t = 20 #Número de periodos en el futuro
df = 100 #fitT[["estimate"]][["df"]]
min.Price = 9.75
#forma de aprovechar lamanera en que R crea las matrices
#Se estiman los parámetros de la distribución Normal.
media <- mean(l1.tv)
des <- sd(l1.tv)
#Se crea una matriz en la que se almacenan los valores simulados de la v.a.
#asumiendo independencia y estacionariedad. Cuando se crea una matriz se declaran tres parámetros 
#1. valores que contendrá la matriz. 
#2. número de filas.
#3. número de columnas.
#En este caso aprovechamos para realizar las simulaciones utilizando la función rnorm() al mismo momento que declaramos la matriz, ya que estableceremos que los valores con los que se debe llenar la matriz son las simulaciones de la v.a.
library(metRology)
## 
## Attaching package: 'metRology'
## The following object is masked from 'package:car':
## 
##     data.ellipse
## The following objects are masked from 'package:base':
## 
##     cbind, rbind
MC.t100 <- matrix(rt.scaled(N*t,df,media,des),t,N)
#se obtiene el valor de los rendimientos acumulados, con lo que se forman las trayectorias de los rendimientos.
MC.t100.Tra <- apply(MC.t100,2,cumsum)
#Se obtiene el último precio de la serie (en este caso el valor 610 de la serie)
P0 <- as.numeric(last(TV1))
#Se genera un vector N con el último precio de la serie.
P0.vector <- rep(P0,N)
#debido a que los rendimientos se obtivieron mediante la diferencia logarítmica 
#manera de reescalar al precio es por medio de la exponencial
MC.t100.Price <- P0*exp(MC.t100.Tra)
#Se agrega en la primera fila de la matriz un vector con P0 para tener la referencia
#de inicio de las simulaciones.
MC.t100.Price <- rbind(P0.vector,MC.t100.Price)
#Se agrega una columna con el número de periodos de tiempo en los que se observará 
#el precio 
MC.t100.Price <- cbind(Time=1:(t+1),MC.t100.Price)
#Se  reornenan los datos para graficarlos
MC.t100.Price.DF <- as.data.frame(MC.t100.Price)
MC.t100.Price.tidy <- MC.t100.Price.DF %>% 
  gather(-Time,key = Sim,value = Price)

MC.t100.Price.tidy %>% 
  ggplot(aes(x=Time,
             y=Price,
             group=Sim))+
  geom_line(alpha=0.1,
            color="blue")+
  theme(legend.position = 'none')+
  labs(title = 'MC t100 Simulation',
       subtitle = 'Televisa/USD',
       x="Days",
       y="Price")+
  geom_hline(yintercept = min.Price,
             col="red",
             size=.2)+
  theme_bw()

# Cálculo de la probabilidad de que el precio de Televisa esté por encima del un umbral
# Se dine el objeto Probabilidad
Prob.t100 <- NULL
#
for (i in 2:t){
  Prob.t100[i] <- sum(MC.t100.Price.DF[i,-1]>=min.Price)/N
}

Prob.t100.DF <- tibble(Dias=1:t,
                       Prob.t100)

Prob.t100.DF %>% 
  ggplot(aes(x=Dias,
             y=Prob.t100))+
  geom_line(col="darkblue",
            size=1)+
  labs(title = "Probability of profit",
       x="Days",
       y="Probability")+
  theme_bw()
## Warning: Removed 1 row(s) containing missing values (geom_path).

#Estimación del Valor Esperado
#Se crea el objeto valor esperado
EV.t100 = NA
#Se estima el valor esperado para cada t
for(i in 1:(t+1)){
  EV.t100[i] <- sum(MC.t100.Price.DF[i,-1])/N
}
#Se crea un Data Frame con los Valores esperados estimados
EV.t100.DF <- tibble(Time = 1:(t+1), Esperado =  EV.t100)
#Se genera nuevamente la gráfica de trayectorias agregando una línea de Valor Esperado
MC.t100.Price.tidy %>% 
  ggplot(aes(x=Time,
             y=Price,
             group=Sim))+
  geom_line(alpha=0.1,
            color="blue")+
  theme(legend.position = 'none')+
  labs(title = 'MC t100 Simulation',
       subtitle = 'Televisa/USD',
       x="Days",
       y="Price")+
  geom_hline(yintercept = min.Price,
             col="red",
             size=.2)+
  geom_line(data = EV.t100.DF,inherit.aes = F, aes(x= Time, y=Esperado), size =1)+
  theme_bw()

#Estimación de quantiles para un nivel de confianza del 95%
MC.t100.Qs = NA
MC.t100.Qi = NA

for(i in 1:(t+1)){
  MC.t100.Qi[i] <- as.numeric(quantile(MC.t100.Price.DF[i,-1], c(0.025,0.975))[1])
  MC.t100.Qs[i] <- as.numeric(quantile(MC.t100.Price.DF[i,-1], c(0.025,0.975))[2])
}
## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames
MC.t100.BandasQ <- tibble(Time =1:(t+1), sup = MC.t100.Qs, inf = MC.t100.Qi)
MC.t100.Price.tidy %>% 
  ggplot(aes(x=Time,
             y=Price,
             group=Sim))+
  geom_line(alpha=0.1,
            color="blue")+
  theme(legend.position = 'none')+
  labs(title = 'MC t100 Simulation',
       subtitle = 'Televisa/USD',
       x="Days",
       y="Price")+
  geom_hline(yintercept = min.Price,
             col="red",
             size=.2)+
  geom_ribbon(data = MC.t100.BandasQ,inherit.aes = F, aes(x= Time, ymin = inf, ymax = sup), fill = "light blue", alpha = 0.6)+
  geom_line(data = EV.t100.DF,inherit.aes = F, aes(x= Time, y=Esperado), size =1)+
  theme_bw()

#Se almacenan los Valores Esparados estimados con la t100
Resultados[3,] <- EV.t100

———————-Metodo Monte Carlo (t-Student)——————

N = 10000 #número de realizaciones, estados de la naturaleza, trayectorias.
t = 20 #Número de periodos en el futuro
#fitT[["estimate"]][["df"]]
min.Price = 9.75 
#forma de aprovechar lamanera en que R crea las matrices
#Se estiman los parámetros de la distribución Normal.
media <- mean(l1.tv)
des <- sd(l1.tv)
#install.packages("fitdistrplus")
#Se estiman los parámetros de la distribución t-student
#Con este método se logra un mejor ajuste ya que esta implementación de la distribución t considera además de los DF la media y la SD como parametros.
library(fitdistrplus)
## Loading required package: survival
library(metRology)
l1.tv<-as.numeric(l1.tv$TV.Adjusted)
t.df <- fitdist(l1.tv, dist="t.scaled",start = list(df=5,mean=mean(l1.tv),sd=sd(l1.tv)))
t5.df <- t.df[["estimate"]][["df"]]
#Se cea una matriz en la que se almacenan los valores simulados de la v.a.
#asumiendo independencia y estacionariedad.
#cuando se crea una matriz se declaran tres parámetros 
#1. valores que contendrá la matriz. 
#2. número de filas.
#3. número de columnas.
#En este caso aprovechamos para realizar las simulaciones utilizando la función 
#rnorm() al mismo momento que declaramos la matriz, ya que estableceremos
#que los valores con los que se debe llenar la matriz son las simulaciones de 
# la v.a.
MC.t5 <- matrix(rt.scaled(N*t,t5.df,media,des),t,N)
#se obtiene el valor de los rendimientos acumulados, con lo que se forman las 
#trayectorias de los rendimientos.
MC.t5.Tra <- apply(MC.t5,2,cumsum)
#Se obtiene el último precio de la serie (en este caso el valor 346 de la serie)
P0 <- as.numeric(last(TV1))
#Se genera un vector N con el último precio de la serie.
P0.vector <- rep(P0,N)
#debido a que los rendimientos se obtivieron mediante la diferencia logarítmica manera de reescalar al precio es por medio de la exponencial
MC.t5.Price <- P0*exp(MC.t5.Tra)
#Se agrega en la primera fila de la matriz un vector con P0 para tener la referencia de inicio de las simulaciones.
MC.t5.Price <- rbind(P0.vector,MC.t5.Price)
#Se agrega una columna con el número de periodos de tiempo en los que se observará el precio 
MC.t5.Price <- cbind(Time=1:(t+1),MC.t5.Price)
max(TV1)
## [1] 14.57111
#Se  reornenan los datos para graficarlos
MC.t5.Price.DF <- as.data.frame(MC.t5.Price)
MC.t5.Price.tidy <- MC.t5.Price.DF %>% 
  gather(-Time,key = Sim,value = Price)

MC.t5.Price.tidy %>% 
  ggplot(aes(x=Time,
             y=Price,
             group=Sim))+
  geom_line(alpha=0.1,
            color="blue")+
  theme(legend.position = 'none')+
  labs(title = 'MC t5 Simulation',
       subtitle = 'TELEVISA/USD',
       x="Days",
       y="Price")+
  geom_hline(yintercept = min.Price,
             col="red",
             size=.2)+
  theme_bw()

# Cálculo de la probabilidad de que el precio de Televisa esté por encima del un umbral. Se dine el objeto Probabilidad
Prob.t5 <- NULL
#
for (i in 2:t){
  Prob.t5[i] <- sum(MC.t5.Price.DF[i,-1]>=min.Price)/N
}

Prob.t5.DF <- tibble(Dias=1:t,
                       Prob.t5)

Prob.t5.DF %>% 
  ggplot(aes(x=Dias,
             y=Prob.t5))+
  geom_line(col="darkblue",
            size=1)+
  labs(title = "Probability of profit",
       x="Days",
       y="Probability")+
  theme_bw()
## Warning: Removed 1 row(s) containing missing values (geom_path).

#Estimación del Valor Esperado
#Se crea el objeto valor esperado
EV.t5 = NA
#Se estima el valor esperado para cada t
for(i in 1:(t+1)){
  EV.t5[i] <- sum(MC.t5.Price.DF[i,-1])/N
}
#Se crea un Data Frame con los Valores esperados estimados
EV.t5.DF <- tibble(Time = 1:(t+1), Esperado =  EV.t5)
#Se genera nuevamente la gráfica de trayectorias agregando una linea de Valor Esperado
MC.t5.Price.tidy %>% 
  ggplot(aes(x=Time,
             y=Price,
             group=Sim))+
  geom_line(alpha=0.1,
            color="blue")+
  theme(legend.position = 'none')+
  labs(title = 'MC t5 Simulation',
       subtitle = 'Televisa/USD',
       x="Days",
       y="Price")+
  geom_hline(yintercept = min.Price,
             col="red",
             size=.2)+
  geom_line(data = EV.t5.DF,inherit.aes = F, aes(x= Time, y=Esperado), size =1)+
  theme_bw()

#Estimación de quantiles para un nivel de confianza del 95%
MC.t5.Qs = NA
MC.t5.Qi = NA

for(i in 1:(t+1)){
  MC.t5.Qi[i] <- as.numeric(quantile(MC.t5.Price.DF[i,-1], c(0.025,0.975))[1])
  MC.t5.Qs[i] <- as.numeric(quantile(MC.t5.Price.DF[i,-1], c(0.025,0.975))[2])
}
## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames
MC.t5.BandasQ <- tibble(Time =1:(t+1), sup = MC.t5.Qs, inf = MC.t5.Qi)
MC.t5.Price.tidy %>% 
  ggplot(aes(x=Time,
             y=Price,
             group=Sim))+
  geom_line(alpha=0.1,
            color="blue")+
  theme(legend.position = 'none')+
  labs(title = 'MC t5 Simulation',
       subtitle = 'TELEVISA/USD',
       x="Days",
       y="Price")+
  geom_hline(yintercept = min.Price,
             col="red",
             size=.2)+
  geom_ribbon(data = MC.t5.BandasQ,inherit.aes = F, aes(x= Time, ymin = inf, ymax = sup), fill = "red", alpha = 0.6)+
  geom_line(data = EV.t5.DF,inherit.aes = F, aes(x= Time, y=Esperado), size =1)+
  theme_bw()

#Se almacenan los Valores Esparados estimados con la t100
Resultados[4,] <- EV.t5

———————–Metodo Monte Carlo (NIG)—————–

N = 10000 #número de realizaciones, estados de la naturaleza, trayectorias.
t = 20 #Número de periodos en el futuro
#fitT[["estimate"]][["df"]]
min.Price = 9.75 
#forma de aprovechar lamanera en que R crea las matrices
#Se estiman los parámetros de la distribución NIG
#Cco este método se logra un mejor ajuste ya que esta implementación de la 
#distribución t considera además de los DF la media y la SD como parametros.
library(GeneralizedHyperbolic)
fit.tv <- nigFit(l1.tv)
fit.tv
## 
## Data:      l1.tv 
## Parameter estimates:
##         mu       delta       alpha        beta  
## -0.0005745   0.0277988  22.6190296   0.0083158  
## Likelihood:         1269.709 
## criterion :         MLE 
## Method:             Nelder-Mead 
## Convergence code:   0 
## Iterations:         213
#install.packages("GeneralizedHyperbolic")
#Se cea una matriz en la que se almacenan los valores simulados de la v.a.
#asumiendo independencia y estacionariedad.
#cuando se crea una matriz se declaran tres parámetros 
#1. valores que contendrá la matriz. 
#2. número de filas.
#3. número de columnas.
#En este caso aprovechamos para realizar las simulaciones utilizando la función 
#rnorm() al mismo momento que declaramos la matriz, ya que estableceremos
#que los valores con los que se debe llenar la matriz son las simulaciones de 
# la v.a.
MC.NIG <- matrix(rnig(length(N*t),mu = fit.tv$param[1], delta = fit.tv$param[2], alpha = fit.tv$param[3], beta = fit.tv$param[4]),t,N)
#Forma paso a paso de hacer el metodo Montecarlo
MC.NIG <- matrix(0,t,N)
for(i in 1:t){
  for(j in 1:N){
    MC.NIG[i,j] <- rnig(1,mu = fit.tv$param[1], delta = fit.tv$param[2], alpha = fit.tv$param[3], beta = fit.tv$param[4])
  }
}
#se obtiene el valor de los rendimientos acumulados, con lo que se forman las 
#trayectorias de los rendimientos.
MC.NIG.Tra <- apply(MC.NIG,2,cumsum)
#Se obtiene el último precio de la serie (en este caso el valor 346 de la serie)
P0 <- as.numeric(last(TV1))
#Se genera un vector N con el último precio de la serie.
P0.vector <- rep(P0,N)
#debido a que los rendimientos se obtivieron mediante la diferencia logarítmica 
#manera de reescalar al precio es por medio de la exponencial
MC.NIG.Price <- P0*exp(MC.NIG.Tra)
#Se agrega en la primera fila de la matriz un vector con P0 para tener la referencia
#de inicio de las simulaciones.
MC.NIG.Price <- rbind(P0.vector,MC.NIG.Price)
#Se agrega una columna con el número de periodos de tiempo en los que se observará 
#el precio 
MC.NIG.Price <- cbind(Time=1:(t+1),MC.NIG.Price)
#Se  reornenan los datos para graficarlos
MC.NIG.Price.DF <- as.data.frame(MC.NIG.Price)
MC.NIG.Price.tidy <- MC.NIG.Price.DF %>% 
  gather(-Time,key = Sim,value = Price)

MC.NIG.Price.tidy %>% 
  ggplot(aes(x=Time,
             y=Price,
             group=Sim))+
  geom_line(alpha=0.1,
            color="blue")+
  theme(legend.position = 'none')+
  labs(title = 'MC NIG Simulation',
       subtitle = 'TELEVISA/USD',
       x="Days",
       y="Price")+
  geom_hline(yintercept = min.Price,
             col="red",
             size=.2)+
  theme_bw()

# Cálculo de la probabilidad de que el precio del BTC esté por encima del un umbral
# Se dine el objeto Probabilidad
Prob.NIG <- NULL
#
for (i in 2:t){
  Prob.NIG[i] <- sum(MC.NIG.Price.DF[i,-1]>=min.Price)/N
}

Prob.NIG.DF <- tibble(Dias=1:t,
                     Prob.NIG)

Prob.NIG.DF %>% 
  ggplot(aes(x=Dias,
             y=Prob.NIG))+
  geom_line(col="darkblue",
            size=1)+
  labs(title = "Probability of profit",
       x="Days",
       y="Probability")+
  theme_bw()
## Warning: Removed 1 row(s) containing missing values (geom_path).

#Estimación del Valor Esperado
#Se crea el objeto valor esperado
EV.NIG = NA
#Se estima el valor esperado para cada t
for(i in 1:(t+1)){
  EV.NIG[i] <- sum(MC.NIG.Price.DF[i,-1])/N
}
#Se crea un Data Frame con los Valores esperados estimados
EV.NIG.DF <- tibble(Time = 1:(t+1), Esperado =  EV.NIG)
#Se genera nuevamente la gráfica de trayectorias agregando una linea de Valor 
#Esperado
MC.NIG.Price.tidy %>% 
  ggplot(aes(x=Time,
             y=Price,
             group=Sim))+
  geom_line(alpha=0.1,
            color="blue")+
  theme(legend.position = 'none')+
  labs(title = 'MC NIG Simulation',
       subtitle = 'TELEVISA/USD',
       x="Days",
       y="Price")+
  geom_hline(yintercept = min.Price,
             col="red",
             size=.2)+
  geom_line(data = EV.NIG.DF,inherit.aes = F, aes(x= Time, y=Esperado), size =1)+
  theme_bw()

#Estimación de quantiles para un nivel de confianza del 95%
MC.NIG.Qs = NA
MC.NIG.Qi = NA

for(i in 1:(t+1)){
  MC.NIG.Qi[i] <- as.numeric(quantile(MC.NIG.Price.DF[i,-1], c(0.025,0.975))[1])
  MC.NIG.Qs[i] <- as.numeric(quantile(MC.NIG.Price.DF[i,-1], c(0.025,0.975))[2])
}
## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames
MC.NIG.BandasQ <- tibble(Time =1:(t+1), sup = MC.NIG.Qs, inf = MC.NIG.Qi)
MC.NIG.Price.tidy %>% 
  ggplot(aes(x=Time,
             y=Price,
             group=Sim))+
  geom_line(alpha=0.1,
            color="blue")+
  theme(legend.position = 'none')+
  labs(title = 'MC NIG Simulation',
       subtitle = 'TELEVISA/USD',
       x="Days",
       y="Price")+
  geom_hline(yintercept = min.Price,
             col="red",
             size=.2)+
  geom_ribbon(data = MC.NIG.BandasQ,inherit.aes = F, aes(x= Time, ymin = inf, ymax = sup), fill = "red", alpha = 0.6)+
  geom_line(data = EV.NIG.DF,inherit.aes = F, aes(x= Time, y=Esperado), size =1)+
  theme_bw()

#Se almacenan los Valores Esparados estimados con la t100
Resultados[5,] <- EV.NIG

——————–Metodo Monte Carlo (Bootstrap)—————–

N = 10000 #número de realizaciones, estados de la naturaleza, trayectorias.
t = 20 #Número de periodos en el futuro
#fitT[["estimate"]][["df"]]
min.Price = 9.75
#forma de aprovechar lamanera en que R crea las matrices
#Se hace el muestreo con remuestreo
Btstrp <- matrix(sample(l1.tv,N*t,replace = T),t,N)
#se obtiene el valor de los rendimientos acumulados, con lo que se forman las 
#trayectorias de los rendimientos.
Btstrp.Tra <- apply(Btstrp,2,cumsum)
#Se obtiene el último precio de la serie (en este caso el valor 610 de la serie)
P0 <- as.numeric(last(TV1))
#Se genera un vector N con el último precio de la serie.
P0.vector <- rep(P0,N)
#debido a que los rendimientos se obtivieron mediante la diferencia logarítmica 
#manera de reescalar al precio es por medio de la exponencial
Btstrp.Price <- P0*exp(Btstrp.Tra)
#Se agrega en la primera fila de la matriz un vector con P0 para tener la referencia
#de inicio de las simulaciones.
Btstrp.Price <- rbind(P0.vector,Btstrp.Price)
#Se agrega una columna con el número de periodos de tiempo en los que se observará 
#el precio 
Btstrp.Price <- cbind(Time=1:(t+1),Btstrp.Price)
#Se  reornenan los datos para graficarlos
Btstrp.Price.DF <- as.data.frame(Btstrp.Price)
Btstrp.Price.tidy <- Btstrp.Price.DF %>% 
  gather(-Time,key = Sim,value = Price)

Btstrp.Price.tidy %>% 
  ggplot(aes(x=Time,
             y=Price,
             group=Sim))+
  geom_line(alpha=0.1,
            color="blue")+
  theme(legend.position = 'none')+
  labs(title = 'Btstrp Simulation',
       subtitle = 'TELEVISA/USD',
       x="Days",
       y="Price")+
  geom_hline(yintercept = min.Price,
             col="red",
             size=.2)+
  theme_bw()

# Cálculo de la probabilidad de que el precio del BTC esté por encima del un umbral
# Se dine el objeto Probabilidad
Prob.Btstrp <- NULL
#
for (i in 2:t){
  Prob.Btstrp[i] <- sum(Btstrp.Price.DF[i,-1]>=min.Price)/N
}

Prob.Btstrp.DF <- tibble(Dias=1:t,
                      Prob.Btstrp)

Prob.Btstrp.DF %>% 
  ggplot(aes(x=Dias,
             y=Prob.Btstrp))+
  geom_line(col="darkblue",
            size=1)+
  labs(title = "Probability of profit",
       x="Days",
       y="Probability")+
  theme_bw()
## Warning: Removed 1 row(s) containing missing values (geom_path).

#Estimación del Valor Esperado
#Se crea el objeto valor esperado
EV.Btstrp = NA
#Se estima el valor esperado para cada t
for(i in 1:(t+1)){
  EV.Btstrp[i] <- sum(Btstrp.Price.DF[i,-1])/N
}
#Se crea un Data Frame con los Valores esperados estimados
EV.Btstrp.DF <- tibble(Time = 1:(t+1), Esperado =  EV.Btstrp)
#Se genera nuevamente la gráfica de trayectorias agregando una linea de Valor 
#Esperado
Btstrp.Price.tidy %>% 
  ggplot(aes(x=Time,
             y=Price,
             group=Sim))+
  geom_line(alpha=0.1,
            color="blue")+
  theme(legend.position = 'none')+
  labs(title = 'Btstrp Simulation',
       subtitle = 'TELEVISA/USD',
       x="Days",
       y="Price")+
  geom_hline(yintercept = min.Price,
             col="red",
             size=.2)+
  geom_line(data = EV.Btstrp.DF,inherit.aes = F, aes(x= Time, y=Esperado), size =1)+
  theme_bw()

#Estimación de quantiles para un nivel de confianza del 95%
Btstrp.Qs = NA
Btstrp.Qi = NA

for(i in 1:(t+1)){
  Btstrp.Qi[i] <- as.numeric(quantile(Btstrp.Price.DF[i,-1], c(0.025,0.975))[1])
  Btstrp.Qs[i] <- as.numeric(quantile(Btstrp.Price.DF[i,-1], c(0.025,0.975))[2])
}
## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames
Btstrp.BandasQ <- tibble(Time =1:(t+1), sup = Btstrp.Qs, inf = Btstrp.Qi)
Btstrp.Price.tidy %>% 
  ggplot(aes(x=Time,
             y=Price,
             group=Sim))+
  geom_line(alpha=0.1,
            color="blue")+
  theme(legend.position = 'none')+
  labs(title = 'Btstrp Simulation',
       subtitle = 'TELEVISA/USD',
       x="Days",
       y="Price")+
  geom_hline(yintercept = min.Price,
             col="red",
             size=.2)+
  geom_ribbon(data = Btstrp.BandasQ,inherit.aes = F, aes(x= Time, ymin = inf, ymax = sup), fill = "red", alpha = 0.6)+
  geom_line(data = EV.Btstrp.DF,inherit.aes = F, aes(x= Time, y=Esperado), size =1)+
  theme_bw()

#Se almacenan los Valores Esparados estimados con la t100
Resultados[6,] <- EV.Btstrp
#install.packages("MLmetrics")
#Estimación de Métricas.
# para determinar que modelo nos permite reliazar un amejor estimación del precio, utilizaremos 
#el MSE, RMSE, MAPE
library(MLmetrics)
## 
## Attaching package: 'MLmetrics'
## The following object is masked from 'package:base':
## 
##     Recall
Metricas <- matrix(0,5,3)
empirica <- c(P0,TV.observed)
Metricas[1,1] <- MSE(EV.Norm, empirica)
Metricas[1,2] <- RMSE(EV.Norm, empirica)
Metricas[1,3] <- MAPE(EV.Norm, empirica)
Metricas[2,1] <- MSE(EV.t100, empirica)
Metricas[2,2] <- RMSE(EV.t100, empirica)
Metricas[2,3] <- MAPE(EV.t100, empirica)
Metricas[3,1] <- MSE(EV.t5, empirica)
Metricas[3,2] <- RMSE(EV.t5, empirica)
Metricas[3,3] <- MAPE(EV.t5, empirica)
Metricas[4,1] <- MSE(EV.NIG, empirica)
Metricas[4,2] <- RMSE(EV.NIG, empirica)
Metricas[4,3] <- MAPE(EV.NIG, empirica)
Metricas[5,1] <- MSE(EV.Btstrp, empirica)
Metricas[5,2] <- RMSE(EV.Btstrp, empirica)
Metricas[5,3] <- MAPE(EV.Btstrp, empirica)
Metricas
##           [,1]      [,2]       [,3]
## [1,] 0.3684198 0.6069759 0.04721366
## [2,] 0.3763009 0.6134337 0.04790645
## [3,] 0.3531900 0.5942979 0.04553086
## [4,] 0.3794010 0.6159553 0.04851950
## [5,] 0.3785945 0.6153003 0.04831935