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)

Ajuste de Distribución y simulaciones Monte Carlo y Bootstrap

1. Gráfica de serie de datos.

  • Serie de datos *
SPWR<-getSymbols("SPWR", from="2017-07-01", to="2022-06-30", 
                    auto.assign = F)[,6]
  • Gráfica *
ggplot(SPWR, aes(x=1:nrow(SPWR), y=SPWR.Adjusted))+
  geom_line(col="darkorchid3",
            size=1)+
  labs(title="Serie de precios SPWR",
       y="Precio",
       x="Observaciones",
       caption = "Fuente Yahoo Finance") +
   theme(plot.title = element_text(color = "deeppink", size = 18, hjust = 0.5),
         axis.title.x.bottom = element_text(color = "deeppink"),
         axis.title.y.left= element_text(color = "deeppink"),
         panel.background = element_rect(fill="white", color = "blue3"))

2. Gráfica de densidad de la serie de datos.

#Se obtiene la densidad de la serie
ggplot(SPWR, aes(x=SPWR.Adjusted))+
  geom_density(color= "red",
               fill = "grey",
               alpha = 0.5)+
  labs(title="Gráfica de densidad",
       y="Density",
       x="SPWR precio",
       caption = "Fuente Yahoo Finance")+
   theme(plot.title = element_text(color = "deeppink", size = 18, hjust = 0.5),
         axis.title.x.bottom = element_text(color = "deeppink"),
         axis.title.y.left= element_text(color = "deeppink"),
         panel.background = element_rect(fill="white", color = "blue3"))

3. Resultado de la prueba estadística ADF.

#Se estima la prueba Augemnted Dickey-Fuller
spwr.adf.test <- adf.test(SPWR)
spwr.adf.test
## 
##  Augmented Dickey-Fuller Test
## 
## data:  SPWR
## Dickey-Fuller = -2.3495, Lag order = 10, p-value = 0.4304
## alternative hypothesis: stationary

Como el p-value es >.05, entonces la serie tiene raíz unitaria, por tanto, no es estacionaria.

4. Mencionar si fue necesario obtener rendimientos logarítmicos de la serie y dar una explicación.

En este caso, es necesario obtener los rendimientos logarítmicos con la finalidad de inducir estacionariedad y que pueda ser modelada. Una ventaja por la que se suelen preferir los rendimientos logarítmicos es la forma de agregación de los periodos.

#En esta linea de código se estiman las diferencias logarítmicas aplicando de manera anidada las funciones diff() y log()
spwr.L1.return <- diff(log(SPWR))
#se almacenan las diferencias logarítminas en un data frame para poder graficarlas mediante ggplot2()
spwr.L1.return.df <- as.data.frame(spwr.L1.return)
ggplot(spwr.L1.return.df, aes(x=1:nrow(spwr.L1.return.df),y=spwr.L1.return))+
  geom_line(col="darkorchid3",
            size=1)+
  labs(title="Retornos Logarítmicos SPWR",
       y="Rendimiento",
       x="Observaciones",
       caption = "Fuente Yahoo Finance")+
   theme(plot.title = element_text(color = "deeppink", size = 18, hjust = 0.5),
         axis.title.x.bottom = element_text(color = "deeppink"),
         axis.title.y.left= element_text(color = "deeppink"),
         panel.background = element_rect(fill="white", color = "blue3"))
## Don't know how to automatically pick scale for object of type xts/zoo. Defaulting to continuous.
## Warning: Removed 1 row(s) containing missing values (geom_path).

** ADF LOGARÍTMICA **

#Se estima nuevamente la prueba Dickey-Fuller para determinar si la diferenciación da estacionariedad a la serie.
spwr.L1.return <- na.omit(spwr.L1.return)

spwr.L1.return.adf.test <- adf.test(spwr.L1.return)
## Warning in adf.test(spwr.L1.return): p-value smaller than printed p-value
spwr.L1.return.adf.test 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  spwr.L1.return
## Dickey-Fuller = -10.752, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary

Como observamos, p-value<.05, por tanto la serie es estacionaria.

#Se obtiene la densidad de la serie
ggplot(spwr.L1.return, aes(x=SPWR.Adjusted))+
  geom_density(color= "red",
               fill = "grey",
               alpha = 0.5)+
  labs(title="Gráfica de densidad",
       y="Density",
       x="SPWR precio",
       caption = "Fuente Yahoo Finance")+
   theme(plot.title = element_text(color = "deeppink", size = 18, hjust = 0.5),
         axis.title.x.bottom = element_text(color = "deeppink"),
         axis.title.y.left= element_text(color = "deeppink"),
         panel.background = element_rect(fill="white", color = "blue3"))

5. Estadística descriptiva de la serie, dando una breve explicación de la media, desviación estándar, sesgo y kurtosis.

#Se estiman la estadística descriptiva de las series de rendimientos.
mean(spwr.L1.return)
## [1] 0.0007165092
median(spwr.L1.return)
## [1] 0
var(spwr.L1.return)
##               SPWR.Adjusted
## SPWR.Adjusted    0.00224832
sd(spwr.L1.return)
## [1] 0.04741645
skewness(spwr.L1.return)
## SPWR.Adjusted 
##    0.05769633
kurtosis(spwr.L1.return)
## SPWR.Adjusted 
##      6.047056
qqnorm(spwr.L1.return$SPWR.Adjusted, pch = 1, col="violet", frame = FALSE)
qqline(spwr.L1.return$SPWR.Adjusted, col = "steelblue", lwd = 2)

#Se obtiene el histograma de la serie
ggplot(spwr.L1.return, aes(x=SPWR.Adjusted))+
  geom_histogram(color= "blue",
               fill = "grey",
               alpha = 0.5)+
  geom_vline(aes(xintercept=mean(spwr.L1.return$SPWR.Adjusted)),
             color="red")+
  geom_vline(aes(xintercept=median(spwr.L1.return$SPWR.Adjusted)),
             color="green")+
    labs(title="Histograma",
       y="Count",
       x="SPWR precio",
       caption = "Fuente Yahoo Finance")+
   theme(plot.title = element_text(color = "deeppink", size = 18, hjust = 0.5),
         axis.title.x.bottom = element_text(color = "deeppink"),
         axis.title.y.left= element_text(color = "deeppink"),
         panel.background = element_rect(fill="white", color = "blue3"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Debido a que la kurtosis es >3, la distribución muestral es apenas leptocúrtica. Existe prácticamente nada de asimetría como vemos más claro en el histograma y parece ser no-normal, kurtosis mayor a 3 y varianza diferente de 0.

6. Resultado de la prueba de normalidad Jarque-Bera y una breve explicación.

jarque.bera.test(spwr.L1.return)
## 
##  Jarque Bera Test
## 
## data:  spwr.L1.return
## X-squared = 486.59, df = 2, p-value < 2.2e-16

Como p-value<.05, significa que la distribución es no-normal y lo habíamos observado anteriormente en la gráfica.

7. Parámetros de la distribución Normal que ajustan los datos.

library("GeneralizedHyperbolic")
#obtención de parámetros de la distribucion empirica
fit.spwr <- nigFit(spwr.L1.return)
#Se realiza el ajuste de la distribución con los parámetros obtenidos
adj.spwr <- rnig(length(spwr.L1.return),mu=fit.spwr$param[1],delta=fit.spwr$param[2],alpha =fit.spwr$param[3],beta =fit.spwr$param[4])
spwr.param.mu <- fit.spwr$param[1]
spwr.param.delta <- fit.spwr$param[2]
spwr.param.alpha <- fit.spwr$param[3]
spwr.param.beta <- fit.spwr$param[4]
spwr.param.mu 
##            mu 
## -0.0006063272
spwr.param.delta 
##      delta 
## 0.04241852
spwr.param.alpha
##    alpha 
## 18.90821
spwr.param.beta 
##      beta 
## 0.5922393

8. Gráfica de las simulaciones generadas, valor esperado y nivel de confianza de la simulación Monte Carlo (Normal) y simulación Bootstrap.

min(SPWR$SPWR.Adjusted)
## [1] 2.940406
max(SPWR$SPWR.Adjusted)
## [1] 54.01
mean(SPWR$SPWR.Adjusted)
## [1] 12.10203
median(SPWR$SPWR.Adjusted)
## [1] 6.11002
#se divide la serie en dos subconjuntos, los primeros 1237 valores se utilizarán para realizar el forecast y los últimos 20 servirán para verificar la precisión de las simulaciones.
SPWR1 <- head(SPWR, n = 1237)
#Se obtienen los últimos 20 valores del la serie
SPWR.observed <- tail(SPWR, n = 20)
#Se obtienen las diferencias lograrítmicas de la series de ajuste.
l1.spwr <- diff(log(SPWR1))
l1.spwr.df <- as.data.frame(l1.spwr)
l2.spwr<-diff(log(SPWR.observed))
l2.spwr.df <- as.data.frame(l2.spwr)
#se estima la prueba ADF para rechazar H0: se tiene raíz unitaria.
l2.spwr<-na.omit(l2.spwr)
l1.spwr<-na.omit(l1.spwr)
spwr.adf <- adf.test(l1.spwr)
## Warning in adf.test(l1.spwr): p-value smaller than printed p-value
spwr.adf
## 
##  Augmented Dickey-Fuller Test
## 
## data:  l1.spwr
## Dickey-Fuller = -10.694, Lag order = 10, 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. ****——————————————————****

N = 10000 #número de realizaciones, estados de la naturaleza, trayectorias.
t = 20 #Número de periodos en el futuro
min.Price = 15
#forma de aprovechar lamanera en que R crea las matrices
#Se estiman los parámetros de la distribución Normal.
media <- mean(l1.spwr)
des <- sd(l1.spwr)
#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(SPWR1))
#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="steelblue")+
  theme(legend.position = 'none')+
  labs(title = 'MC Normal Simulation',
       subtitle = 'SPWR/USD',
       x="Days",
       y="Price")+
  geom_hline(yintercept = min.Price,
             col="red",
             size=.5)+ 
  theme(plot.title = element_text(color = "deeppink", size = 18, hjust = 0.5),
        plot.subtitle = element_text(color = "deeppink", size = 14, hjust = 0.5),
         axis.title.x.bottom = element_text(color = "deeppink"),
         axis.title.y.left= element_text(color = "deeppink"),
         panel.background = element_rect(fill="white", color = "blue3"))

#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.5,
            color="steelblue")+
  theme(legend.position = 'none')+
  labs(title = 'MC Normal Simulation',
       subtitle = 'SPWR/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), col="black", size =1)+
   theme(plot.title = element_text(color = "deeppink", size = 18, hjust = 0.5),
        plot.subtitle = element_text(color = "deeppink", size = 14, hjust = 0.5),
         axis.title.x.bottom = element_text(color = "deeppink"),
         axis.title.y.left= element_text(color = "deeppink"),
         panel.background = element_rect(fill="white", color = "blue3"))

#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="steelblue")+
  theme(legend.position = 'none')+
  labs(title = 'MC Normal Simulation',
       subtitle = 'SPWR/USD',
       x="Days",
       y="Price")+
  geom_hline(yintercept = min.Price,
             col="red",
             size=.5)+
  geom_ribbon(data = MC.Norm.BandasQ,inherit.aes = F, aes(x= Time, ymin = inf, ymax = sup), fill = "pink", alpha = 0.5)+
  geom_line(data = EV.Norm.DF,inherit.aes = F, aes(x= Time, y=Esperado), size =1)+
  theme(plot.title = element_text(color = "deeppink", size = 18, hjust = 0.5),
        plot.subtitle = element_text(color = "deeppink", size = 14, hjust = 0.5),
         axis.title.x.bottom = element_text(color = "deeppink"),
         axis.title.y.left= element_text(color = "deeppink"),
         panel.background = element_rect(fill="white", color = "blue3"))

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

Método 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 = 15
#forma de aprovechar lamanera en que R crea las matrices
#Se hace el muestreo con remuestreo
Btstrp <- matrix(sample(l1.spwr,N*t,replace = T),t,N)

#Btstrp <-matrix(0,t,N)
#for(i in 1:t){
 # for(j in 1:N){
  #  l2.tv.df[i,j] <- sample(l1.tv,1,replace = T)
#  }
#}
#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 1237 de la serie)
P0 <- as.numeric(last(SPWR1))
#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="steelblue")+
  theme(legend.position = 'none')+
  labs(title = 'Btstrp Simulation',
       subtitle = 'SPWR/USD',
       x="Days",
       y="Price")+
  geom_hline(yintercept = min.Price,
             col="red",
             size=.5)+
  theme(plot.title = element_text(color = "deeppink", size = 18, hjust = 0.5),
        plot.subtitle = element_text(color = "deeppink", size = 14, hjust = 0.5),
         axis.title.x.bottom = element_text(color = "deeppink"),
         axis.title.y.left= element_text(color = "deeppink"),
         panel.background = element_rect(fill="white", color = "blue3"))

# 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(plot.title = element_text(color = "deeppink", size = 18, hjust = 0.5),
        plot.subtitle = element_text(color = "deeppink", size = 14, hjust = 0.5),
         axis.title.x.bottom = element_text(color = "deeppink"),
         axis.title.y.left= element_text(color = "deeppink"),
         panel.background = element_rect(fill="white", color = "blue3"))
## 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="steelblue")+
  theme(legend.position = 'none')+
  labs(title = 'Btstrp Simulation',
       subtitle = 'SPWR/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(plot.title = element_text(color = "deeppink", size = 18, hjust = 0.5),
        plot.subtitle = element_text(color = "deeppink", size = 14, hjust = 0.5),
         axis.title.x.bottom = element_text(color = "deeppink"),
         axis.title.y.left= element_text(color = "deeppink"),
         panel.background = element_rect(fill="white", color = "blue3"))

#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="steelblue")+
  theme(legend.position = 'none')+
  labs(title = 'Btstrp Simulation',
       subtitle = 'SPWR/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 = "purple", alpha = 0.7)+
  geom_line(data = EV.Btstrp.DF,inherit.aes = F, aes(x= Time, y=Esperado), size =1)+
  theme(plot.title = element_text(color = "deeppink", size = 18, hjust = 0.5),
        plot.subtitle = element_text(color = "deeppink", size = 14, hjust = 0.5),
         axis.title.x.bottom = element_text(color = "deeppink"),
         axis.title.y.left= element_text(color = "deeppink"),
         panel.background = element_rect(fill="white", color = "blue3"))

#Se almacenan los Valores Esparados estimados con la t100
Resultados[3,] <- EV.Btstrp
Resultados
##       [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
## [1,] 17.67 17.62000 18.70000 18.38000 18.88000 19.36000 19.04000 18.86000
## [2,] 17.67 17.70967 17.74147 17.77803 17.83220 17.87049 17.90847 17.94576
## [3,] 17.67 17.70432 17.77957 17.89507 18.05101 18.25308 18.50068 18.79796
##          [,9]    [,10]    [,11]    [,12]    [,13]    [,14]    [,15]    [,16]
## [1,] 18.04000 16.71000 16.48000 17.39000 16.07000 17.93000 17.29000 16.65000
## [2,] 17.97918 18.01346 18.03810 18.07939 18.11746 18.16030 18.19515 18.22667
## [3,] 19.14799 19.55864 20.03315 20.58088 21.20705 21.92491 22.74854 23.68949
##         [,17]    [,18]    [,19]    [,20]    [,21]
## [1,] 17.49000 17.36000 18.13000 16.97000 15.14000
## [2,] 18.27318 18.32297 18.35855 18.39208 18.44213
## [3,] 24.77129 26.01252 27.43474 29.07003 30.94210
#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,2,3)
empirica <- c(P0,SPWR.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.Btstrp, empirica)
Metricas[2,2] <- RMSE(EV.Btstrp, empirica)
Metricas[2,3] <- MAPE(EV.Btstrp, empirica)
Metricas
##           [,1]     [,2]       [,3]
## [1,]  1.590341 1.261087 0.05938398
## [2,] 36.558244 6.046341 0.25284551
dimnames(Metricas)<-list(c("MC Normal", "Bootstrap"), c("MSE", "RMSE", "MAPE"))
Metricas
##                 MSE     RMSE       MAPE
## MC Normal  1.590341 1.261087 0.05938398
## Bootstrap 36.558244 6.046341 0.25284551