Untitled
Lectura de la base de datos
library(readxl)
library(tseries)
library(tidyverse)
library(mice)
library(reshape)
library(lubridate)
library(boot)Humedad <- read_excel("Humedad.xlsx",skip = 0,col_names = TRUE)
Humedad$Cumanda <- as.numeric(Humedad$Cumanda)
Humedad$Espoch <- as.numeric(Humedad$Espoch)
Humedad$Multitud <- as.numeric(Humedad$Multitud)
Humedad$Tunshi <- as.numeric(Humedad$Tunshi)Imputación de los datos paquete mice, método pmm
# Imputacion previa de la base de datos
imputed_data_pmm <- mice(Humedad, m=5, maxit=10
, method = "pmm",seed=1996)##
## iter imp variable
## 1 1 Atillo Cumanda Espoch Multitud Tunshi
## 1 2 Atillo Cumanda Espoch Multitud Tunshi
## 1 3 Atillo Cumanda Espoch Multitud Tunshi
## 1 4 Atillo Cumanda Espoch Multitud Tunshi
## 1 5 Atillo Cumanda Espoch Multitud Tunshi
## 2 1 Atillo Cumanda Espoch Multitud Tunshi
## 2 2 Atillo Cumanda Espoch Multitud Tunshi
## 2 3 Atillo Cumanda Espoch Multitud Tunshi
## 2 4 Atillo Cumanda Espoch Multitud Tunshi
## 2 5 Atillo Cumanda Espoch Multitud Tunshi
## 3 1 Atillo Cumanda Espoch Multitud Tunshi
## 3 2 Atillo Cumanda Espoch Multitud Tunshi
## 3 3 Atillo Cumanda Espoch Multitud Tunshi
## 3 4 Atillo Cumanda Espoch Multitud Tunshi
## 3 5 Atillo Cumanda Espoch Multitud Tunshi
## 4 1 Atillo Cumanda Espoch Multitud Tunshi
## 4 2 Atillo Cumanda Espoch Multitud Tunshi
## 4 3 Atillo Cumanda Espoch Multitud Tunshi
## 4 4 Atillo Cumanda Espoch Multitud Tunshi
## 4 5 Atillo Cumanda Espoch Multitud Tunshi
## 5 1 Atillo Cumanda Espoch Multitud Tunshi
## 5 2 Atillo Cumanda Espoch Multitud Tunshi
## 5 3 Atillo Cumanda Espoch Multitud Tunshi
## 5 4 Atillo Cumanda Espoch Multitud Tunshi
## 5 5 Atillo Cumanda Espoch Multitud Tunshi
## 6 1 Atillo Cumanda Espoch Multitud Tunshi
## 6 2 Atillo Cumanda Espoch Multitud Tunshi
## 6 3 Atillo Cumanda Espoch Multitud Tunshi
## 6 4 Atillo Cumanda Espoch Multitud Tunshi
## 6 5 Atillo Cumanda Espoch Multitud Tunshi
## 7 1 Atillo Cumanda Espoch Multitud Tunshi
## 7 2 Atillo Cumanda Espoch Multitud Tunshi
## 7 3 Atillo Cumanda Espoch Multitud Tunshi
## 7 4 Atillo Cumanda Espoch Multitud Tunshi
## 7 5 Atillo Cumanda Espoch Multitud Tunshi
## 8 1 Atillo Cumanda Espoch Multitud Tunshi
## 8 2 Atillo Cumanda Espoch Multitud Tunshi
## 8 3 Atillo Cumanda Espoch Multitud Tunshi
## 8 4 Atillo Cumanda Espoch Multitud Tunshi
## 8 5 Atillo Cumanda Espoch Multitud Tunshi
## 9 1 Atillo Cumanda Espoch Multitud Tunshi
## 9 2 Atillo Cumanda Espoch Multitud Tunshi
## 9 3 Atillo Cumanda Espoch Multitud Tunshi
## 9 4 Atillo Cumanda Espoch Multitud Tunshi
## 9 5 Atillo Cumanda Espoch Multitud Tunshi
## 10 1 Atillo Cumanda Espoch Multitud Tunshi
## 10 2 Atillo Cumanda Espoch Multitud Tunshi
## 10 3 Atillo Cumanda Espoch Multitud Tunshi
## 10 4 Atillo Cumanda Espoch Multitud Tunshi
## 10 5 Atillo Cumanda Espoch Multitud TunshiData_Impu_pmm <- mice::complete(imputed_data_pmm)
Humedad_completa <- Data_Impu_pmmGráficos de las series
Numero <- c(1:1095)
df <- cbind(Numero,Humedad)
df <- melt(df,id.vars = "Numero")
ggplot(df, aes(x = Numero, y = value, color = variable)) +
geom_line()+labs(x="Observaciones",y="Humedad Relativa") +
ggtitle("Series metorológicas") +
theme(plot.title = element_text(hjust = 0.5))
# Histograma de todos los datos de la Humedad Relativa
ggplot(df,aes(x=value)) +
geom_histogram(col="black",fill="pink3",bins = 20) +
labs(x="Humedad Relativa",y="Frecuencia") +
ggtitle("Histograma de los datos de la Humedad Relativa")+
theme(plot.title = element_text(hjust = 0.5))
Transformación de la base de datos a series de tiempo
Numero <- c(1:1095)
df <- cbind(Numero,Humedad_completa)
df <- melt(df,id.vars = "Numero")
d <- ts(df$value)bootf <- function(tsb){ #statistic for time - series bootstrap
fit <- arima(tsb, order = c(1,0,0)) #fi t with AR(1)
return(coef(fit)[1])}
boot2 <- tsboot(d, bootf, R = 2000, l = 16, sim ="fixed" )
mean(boot2$t[,1]) #mean estimate [1] 0.07234437 ## [1] 0.7791205quantile(boot2$t[,1], probs = c(0.025,0.975))## 2.5% 97.5%
## 0.7453390 0.8137052hist(boot2$t[,1], bins=50)## Warning in plot.window(xlim, ylim, "", ...): "bins" is not a graphical parameter## Warning in title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...): "bins"
## is not a graphical parameter## Warning in axis(1, ...): "bins" is not a graphical parameter## Warning in axis(2, ...): "bins" is not a graphical parameter
summary(boot2$t[,1])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.7222 0.7666 0.7795 0.7791 0.7910 0.8373boot.ci(boot.out=boot2,type=c("norm","basic","perc"))## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot2, type = c("norm", "basic", "perc"))
##
## Intervals :
## Level Normal Basic Percentile
## 95% ( 0.8531, 0.9219 ) ( 0.8529, 0.9213 ) ( 0.7453, 0.8138 )
## Calculations and Intervals on Original ScaleView(boot2$t)
View(boot2$t[,1])
print(boot2)##
## BLOCK BOOTSTRAP FOR TIME SERIES
##
## Fixed Block Length of 16
##
## Call:
## tsboot(tseries = d, statistic = bootf, R = 2000, l = 16, sim = "fixed")
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.8333245 -0.05420392 0.01754601base <- boot2$t[,1]
hist(base, freq = F,
ylab = "Densidad",
xlab = "Estadísticos bootstrap", main = "Histogramas de las muestras bootstrap")
dz <- density(base)
lines(dz, col = "red", lwd = 1)
curve(dnorm(x, mean(base), sd(base)),
col = "blue", lwd = 1, add = TRUE)
Bootstrap estacionario (longitud del bloque \(l\sim\)Geometrica)
boot3 <- tsboot(d, bootf, R = 1000, l = 16, sim ="geom" )
mean(boot3$t[,1])## [1] 0.7797182quantile(boot3$t[,1], probs = c(0.025,0.975))## 2.5% 97.5%
## 0.7403281 0.8215009hist(boot3$t[,1])
summary(boot3$t[,1])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.7140 0.7639 0.7791 0.7797 0.7951 0.8529boot.ci(boot.out=boot3,type=c("norm","basic","perc","bca"))## Warning in boot.ci(boot.out = boot3, type = c("norm", "basic", "perc", "bca")):
## BCa intervals not defined for time series bootstraps## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot3, type = c("norm", "basic", "perc", "bca"))
##
## Intervals :
## Level Normal Basic Percentile
## 95% ( 0.8448, 0.9291 ) ( 0.8451, 0.9264 ) ( 0.7402, 0.8215 )
## Calculations and Intervals on Original Scaleprint(boot3)##
## STATIONARY BOOTSTRAP FOR TIME SERIES
##
## Average Block Length of 16
##
## Call:
## tsboot(tseries = d, statistic = bootf, R = 1000, l = 16, sim = "geom")
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.8333245 -0.05360628 0.02149265Normalizar la base Horiginal
HumedadR <- (d-mean(d))/sd(d)Anomalías normalizadas
Anomalias <- HumedadR - boot3$t[,1]## Warning in `-.default`(HumedadR, boot3$t[, 1]): longitud de objeto mayor no es
## múltiplo de la longitud de uno menorhist(Anomalias)
boxplot(Anomalias)
summary(Anomalias)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -8.8529 -1.3522 -0.6252 -0.7797 -0.0944 0.8272k <- length(d)
d[-1] <- d[-1] + d[-k] # introduce auto-correlation
MMB <- function(x, R=100, Tamaniobloque=20){
# Inicialización
n <- length(x)
lagCorMMB <- ct <- numeric(R)
seriesB1 <- numeric(n)
# Función para Bootstrap en bloques moviles
corT <- function(x=tseries,N=n,bl=Tamaniobloque){
for (i in 1:ceiling(N/bl)){
endpoint <- sample(bl:N, size=1)
seriesB1[(i-1)*bl+1:bl] <- x[endpoint-(bl:1)+1]
}
seriesB1 <- seriesB1[1:N]
a <- cor(seriesB1[-1],seriesB1[-N])
return(a)
}
ct <- replicate(R, corT(x))
return(ct)
}
MB <- MMB(d,100,20)
MB[1:10]
mean(MB)
plot(ts(MB))
hist(MB)
quantile(MB, c(0.025,0.975))Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.