R Markdown
# Limpando Environment do R
rm(list = ls())
# Carregando pacotes via Pacman
library(pacman)
pacman::p_load(
"tidyverse",
"readxl",
"lubridate",
"stringr",
"janitor",
"arsenal",
"diffdf",
"knitr",
"kableExtra",
"RQuantLib",
"ggpubr",
"splines",
"stats",
"VGAM",
"stats4")
#Número de Repetições
Rp<-10000
#Tamanho da amostra padrão
sz<- 1000
#Número do grupo
K <- 9
#Amostras
n1<- K*1000
n2<- K*100
n3<- K*10
n4<- K+5
#Matrizes de resultado
TT<-matrix(0,ncol=12,nrow=Rp) # estimates matrix
TS<-matrix(0,ncol=12,nrow=Rp) # Std Err. matrix
############################################################################## QUESTÃO 3
for(i in 1:Rp) {
Y<-rexp(sz,rate=2) # Exponencial
X<-exp(Y) # Pareto com par=2
#Método de Momentos
#MM estimates n1
TT[i,1]<-mean(X)/(mean(X)-1)
TS[i,1]<-( ((TT[i,1]+1)^2 *TT[i,1])/(TT[i,1]+2) )/sz
#MM estimates n2
TT[i,2]<-mean(X[c(1:n2)])/(mean(X[c(1:n2)])-1)
TS[i,2]<-( ((TT[i,2]+1)^2 *TT[i,2])/(TT[i,2]+2) )/n2
#MM estimates n3
TT[i,3]<-mean(X[c(1:n3)])/(mean(X[c(1:n3)])-1)
TS[i,3]<-( ((TT[i,3]+1)^2 *TT[i,3])/(TT[i,3]+2) )/n3
#MM estimates n4
TT[i,4]<-mean(X[c(1:n4)])/(mean(X[c(1:n4)])-1)
TS[i,4]<-( ((TT[i,4]+1)^2 *TT[i,4])/(TT[i,4]+2) )/n4
#Máxima Verossimilhança
#ML estimates n1
TT[i,5]<-1/mean(Y)
TS[i,5]<-TT[i,5]^2/sz
#ML estimates n2
TT[i,6]<-1/mean(Y[c(1:n2)])
TS[i,6]<-TT[i,6]^2/(n2)
#ML estimates n3
TT[i,7]<-1/mean(Y[c(1:n3)])
TS[i,7]<-TT[i,7]^2/(n3)
#ML estimates n4
TT[i,8]<-1/mean(Y[c(1:n4)])
TS[i,8]<-TT[i,8]^2/(n4)
#Mediana
#ME estimates n1
TT[i,9]<-log(2)/log(median(X))
TS[i,9]<-(1/log(median(X))^2) * var(log(X))
#ME estimates n2
TT[i,10]<-log(2)/log(median(X[c(1:n2)]))
TS[i,10]<-(1/log(median(X[c(1:n2)]))^2) * var(log(X[c(1:n2)]))
#ME estimates n3
TT[i,11]<-log(2)/log(median(X[c(1:n3)]))
TS[i,11]<-(1/log(median(X[c(1:n3)]))^2) * var(log(X[c(1:n3)]))
#ME estimates n4
TT[i,12]<-log(2)/log(median(X[c(1:n4)]))
TS[i,12]<-(1/log(median(X[c(1:n4)]))^2) * var(log(X[c(1:n4)]))
}
############################################################################## QUESTÃO 4
############## I)
dTT1 <- density(TT[,1])
dTT2 <- density(TT[,2])
dTT3 <- density(TT[,3])
dTT4 <- density(TT[,4])
dTT5 <- density(TT[,5])
dTT6 <- density(TT[,6])
dTT7 <- density(TT[,7])
dTT8 <- density(TT[,8])
dTT9 <- density(TT[,9])
dTT10 <- density(TT[,10])
dTT11 <- density(TT[,11])
dTT12 <- density(TT[,12])
#Plotando os Gráficos - amostra n5
plot(dTT1, xlim = c(1,3), lty = 1, lwd = 2, main = " ", xlab = " ", ylab = " ")
lines(dTT5, lty = 5, lwd = 2)
lines(dTT9, lty = 9, lwd = 2)
legend("topright", legend = c("MM", "MV", "ME"), lty = 1:3, lwd = 2)
#Plotando os Gráficos - amostra n10
plot(dTT2, xlim = c(1,3), lty = 2, lwd = 2, main = " ", xlab = " ", ylab = " ")
lines(dTT6, lty = 6, lwd = 2)
lines(dTT10, lty = 10, lwd = 2)
legend("topright", legend = c("MM", "MV", "ME"), lty = 1:3, lwd = 2)
#Plotando os Gráficos - amostra n100
plot(dTT3, xlim = c(1,3), lty = 3, lwd = 2, main = " ", xlab = " ", ylab = " ")
lines(dTT7, lty = 7, lwd = 2)
lines(dTT11, lty = 11, lwd = 2)
legend("topright", legend = c("MM", "MV", "ME"), lty = 1:3, lwd = 2)
#Plotando os Gráficos - amostra n1000
plot(dTT4, xlim = c(1,3), lty = 4, lwd = 2, main = " ", xlab = " ", ylab = " ")
lines(dTT8, lty = 8, lwd = 2)
lines(dTT12, lty = 12, lwd = 2)
legend("topright", legend = c("MM", "MV", "ME"), lty = 1:3, lwd = 2)
########### II)
# Convergência em distribuição
TTS <- (TT-2)/sqrt(TS)
summary(TTS[,c(1:4)]) # MM
## V1 V2 V3 V4
## Min. :-18.46929 Min. :-18.12770 Min. :-7.8489 Min. :-2.8543
## 1st Qu.: -0.63236 1st Qu.: -0.61743 1st Qu.:-0.4445 1st Qu.:-0.2297
## Median : 0.21773 Median : 0.23524 Median : 0.2813 Median : 0.3304
## Mean : 0.08991 Mean : 0.09082 Mean : 0.1754 Mean : 0.2817
## 3rd Qu.: 1.01063 3rd Qu.: 1.00478 3rd Qu.: 0.9236 3rd Qu.: 0.8480
## Max. : 4.36681 Max. : 4.35642 Max. : 3.3723 Max. : 2.9597
summary(TTS[,c(5:8)]) # MV
## V1 V2 V3 V4
## Min. :-3.972600 Min. :-3.61727 Min. :-4.10769 Min. :-5.60595
## 1st Qu.:-0.673425 1st Qu.:-0.65400 1st Qu.:-0.65131 1st Qu.:-0.58822
## Median : 0.016242 Median : 0.02099 Median : 0.04322 Median : 0.09831
## Mean : 0.006403 Mean : 0.01014 Mean : 0.00082 Mean : 0.01296
## 3rd Qu.: 0.699443 3rd Qu.: 0.69810 3rd Qu.: 0.68846 3rd Qu.: 0.71820
## Max. : 3.903969 Max. : 4.11300 Max. : 3.24850 Max. : 2.95269
summary(TTS[,c(9:12)]) # ME
## V1 V2 V3
## Min. :-0.3176734 Min. :-0.2607077 Min. :-1.017885
## 1st Qu.:-0.0421525 1st Qu.:-0.0441586 1st Qu.:-0.144964
## Median : 0.0007262 Median : 0.0010536 Median : 0.001520
## Mean : 0.0001809 Mean : 0.0006391 Mean :-0.004906
## 3rd Qu.: 0.0436505 3rd Qu.: 0.0466674 3rd Qu.: 0.142424
## Max. : 0.3028384 Max. : 0.3080538 Max. : 0.803979
## V4
## Min. :-2.80184
## 1st Qu.:-0.40917
## Median :-0.01157
## Mean :-0.02275
## 3rd Qu.: 0.36053
## Max. : 6.50633
# Calculando as densidades do estimador padronizado
densidades_TTS <- list()
for (i in 1:12) {
densidades_TTS[[i]] <- density(TTS[,i])
}
# Plotando as densidades do estimador padronizado
cores <- c("red", "blue", "green")
plot(densidades_TTS[[1]], xlim = range(TTS), lty = 1, lwd = 2, col = cores[1],
main = "Densidades do estimador padronizado", xlab = "TTS", ylab = "Densidade")
for (i in 2:12) {
lines(densidades_TTS[[i]], lty = i%%4, lwd = 2, col = cores[(i-1) %/% 4 + 1])
}
legend("topright", legend = c("MM", "MV", "ME"), lty = 1:3, lwd = 2, col = cores)
############################################################################## QUESTÃO 5
summary(TT[,1])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.208 1.958 2.015 2.010 2.070 2.338
summary(TT[,2])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.191 1.957 2.017 2.011 2.073 2.359
summary(TT[,3])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.025 1.905 2.065 2.068 2.228 3.156
summary(TT[,4])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.069 1.877 2.205 2.308 2.617 9.731
summary(TT[,5])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.777 1.958 2.001 2.002 2.045 2.282
summary(TT[,6])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.785 1.957 2.001 2.003 2.048 2.318
summary(TT[,7])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.396 1.872 2.009 2.023 2.156 3.041
summary(TT[,8])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.8006 1.7283 2.0540 2.1611 2.4751 9.4850
summary(TT[,9])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.627 1.940 2.001 2.004 2.064 2.463
summary(TT[,10])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.677 1.938 2.002 2.005 2.069 2.463
summary(TT[,11])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.109 1.812 2.002 2.031 2.218 3.648
summary(TT[,12])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.6274 1.5623 1.9847 2.1841 2.5838 10.5501
BiasTT<-cbind(mean((TT[,1]-2)), mean((TT[,2]-2)),mean((TT[,3]-2)), mean((TT[,4]-2)),mean((TT[,5]-2)),mean((TT[,6]-2)),mean((TT[,7]-2)),mean((TT[,8]-2)),mean((TT[,9]-2)),mean((TT[,10]-2)),mean((TT[,11]-2)),mean((TT[,12]-2)))
EQMTT<-cbind(mean((TT[,1]-2)^2), mean((TT[,2]-2)^2),mean((TT[,3]-2)^2), mean((TT[,4]-2)^2),mean((TT[,5]-2)^2),mean((TT[,6]-2)^2),mean((TT[,7]-2)^2),mean((TT[,8]-2)^2),mean((TT[,9]-2)^2),mean((TT[,10]-2)^2),mean((TT[,11]-2)^2),mean((TT[,12]-2)^2))
BiasTT
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.009943539 0.01071925 0.06752248 0.30812 0.002396527 0.002875251
## [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 0.02261451 0.1611397 0.00384289 0.004841685 0.03086677 0.1841051
EQMTT
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.008289524 0.009096123 0.06507942 0.4870498 0.004001679 0.00442029
## [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 0.04679186 0.4179535 0.00859376 0.009426885 0.09790079 0.8457834
############################################################################## QUESTÃO 6
summary(TS[,1])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.001837 0.004330 0.004561 0.004549 0.004793 0.006007
summary(TS[,2])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.001991 0.004806 0.005077 0.005059 0.005342 0.006783
summary(TS[,3])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.01544 0.04574 0.05301 0.05380 0.06102 0.11745
summary(TS[,4])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1065 0.2863 0.3847 0.4462 0.5298 6.8234
summary(TS[,5])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.003157 0.003835 0.004004 0.004014 0.004183 0.005206
summary(TS[,6])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.003539 0.004257 0.004451 0.004462 0.004659 0.005969
summary(TS[,7])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.02164 0.03892 0.04485 0.04597 0.05167 0.10278
summary(TS[,8])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.04578 0.21336 0.30134 0.36161 0.43757 6.42606
summary(TS[,9])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.368 1.921 2.076 2.092 2.246 3.428
summary(TS[,10])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.287 1.912 2.075 2.094 2.255 3.338
summary(TS[,11])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.4904 1.5500 1.9981 2.1679 2.6064 16.5608
summary(TS[,12])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.08958 0.90354 1.59560 2.59902 2.92401 97.02166
BiasTS <- cbind(mean((TS[,1]-2)), mean((TS[,2]-2)),mean((TS[,3]-2)), mean((TS[,4]-2)),mean((TS[,5]-2)),mean((TS[,6]-2)),mean((TS[,7]-2)),mean((TS[,8]-2)),mean((TS[,9]-2)),mean((TS[,10]-2)),mean((TS[,11]-2)),mean((TS[,12]-2)))
EQMTS <- cbind(mean((TS[,1]-2)^2), mean((TS[,2]-2)^2),mean((TS[,3]-2)^2), mean((TS[,4]-2)^2),mean((TS[,5]-2)^2),mean((TS[,6]-2)^2),mean((TS[,7]-2)^2),mean((TS[,8]-2)^2),mean((TS[,9]-2)^2),mean((TS[,10]-2)^2),mean((TS[,11]-2)^2),mean((TS[,12]-2)^2))
BiasTS
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] -1.995451 -1.994941 -1.946204 -1.553826 -1.995986 -1.995538 -1.954031
## [,8] [,9] [,10] [,11] [,12]
## [1,] -1.638392 0.09181627 0.09435732 0.1678821 0.5990219
EQMTS
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 3.981825 3.97979 3.787847 2.479572 3.983962 3.982171 3.818334 2.741489
## [,9] [,10] [,11] [,12]
## [1,] 0.06895349 0.07569763 0.8102199 12.74676
############################################################################## QUESTÃO 7
setwd("C:/Users/MarcelAlmeida/Cofco International/FSUSA - Coffee/Marcel/R/")
data_pop <- readxl::read_excel("pop2000.xlsx")
#Características
K = 9
samplesz <- 50000+K*50
samplesz
## [1] 50450
df <- subset(data_pop, pop_2000>samplesz)
df_est <- matrix(0, ncol=3, nrow = 3) # estimates matrix
df_est2 <- matrix(0, ncol=3, nrow = 3) # Std Err. matrix