R Markdown
######################### TRABALHO 1 - ECONOMETRIA 1
######################### PPGE-UFRJ
######################### GRUPO 9 - EMMANUEL & MARCEL
# 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. :-23.98088 Min. :-23.04088 Min. :-7.1276 Min. :-3.1638
## 1st Qu.: -0.67122 1st Qu.: -0.64460 1st Qu.:-0.4209 1st Qu.:-0.2632
## Median : 0.23128 Median : 0.23069 Median : 0.2733 Median : 0.3139
## Mean : 0.05926 Mean : 0.06459 Mean : 0.1719 Mean : 0.2599
## 3rd Qu.: 1.00782 3rd Qu.: 1.01128 3rd Qu.: 0.9070 3rd Qu.: 0.8513
## Max. : 4.06777 Max. : 4.28248 Max. : 3.1926 Max. : 2.6101
summary(TTS[,c(5:8)]) # MV
## V1 V2 V3 V4
## Min. :-4.231820 Min. :-3.715285 Min. :-3.842273 Min. :-4.64995
## 1st Qu.:-0.671731 1st Qu.:-0.665391 1st Qu.:-0.658446 1st Qu.:-0.64487
## Median : 0.013498 Median : 0.012754 Median : 0.037521 Median : 0.08314
## Mean : 0.001488 Mean : 0.004248 Mean : 0.002912 Mean :-0.01209
## 3rd Qu.: 0.682162 3rd Qu.: 0.683896 3rd Qu.: 0.693440 3rd Qu.: 0.71245
## Max. : 3.595819 Max. : 3.624316 Max. : 3.303770 Max. : 2.56738
summary(TTS[,c(9:12)]) # ME
## V1 V2 V3
## Min. :-0.2717270 Min. :-0.2660962 Min. :-1.0125684
## 1st Qu.:-0.0418996 1st Qu.:-0.0440886 1st Qu.:-0.1412299
## Median : 0.0005792 Median : 0.0007247 Median : 0.0022412
## Mean : 0.0001940 Mean : 0.0003200 Mean : 0.0002713
## 3rd Qu.: 0.0428092 3rd Qu.: 0.0447168 3rd Qu.: 0.1473844
## Max. : 0.2278944 Max. : 0.2465801 Max. : 0.7706602
## V4
## Min. :-2.72990
## 1st Qu.:-0.40829
## Median :-0.01479
## Mean :-0.02984
## 3rd Qu.: 0.35623
## Max. : 4.27257
# 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.072 1.956 2.016 2.008 2.070 2.312
summary(TT[,2])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.065 1.955 2.016 2.009 2.074 2.352
summary(TT[,3])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.077 1.910 2.063 2.067 2.224 3.064
summary(TT[,4])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.013 1.860 2.194 2.294 2.620 6.745
summary(TT[,5])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.764 1.958 2.001 2.002 2.044 2.257
summary(TT[,6])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.780 1.957 2.001 2.002 2.047 2.275
summary(TT[,7])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.423 1.870 2.008 2.023 2.158 3.069
summary(TT[,8])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.8918 1.7060 2.0455 2.1470 2.4704 6.3727
summary(TT[,9])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.649 1.941 2.001 2.004 2.063 2.371
summary(TT[,10])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.654 1.938 2.001 2.004 2.066 2.422
summary(TT[,11])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.200 1.816 2.003 2.038 2.226 3.713
summary(TT[,12])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.6201 1.5556 1.9811 2.1757 2.5716 9.3136
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.008345451 0.009331507 0.0670083 0.2937466 0.002083015 0.002475444
## [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 0.02311746 0.1470445 0.003641065 0.004163761 0.03780114 0.1756861
EQMTT
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.009048484 0.009869497 0.0653915 0.4680585 0.003995905 0.004408015
## [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 0.04715633 0.4033023 0.008196453 0.009058549 0.09987027 0.8401895
############################################################################## QUESTÃO 6
summary(TS[,1])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.001498 0.004320 0.004565 0.004543 0.004793 0.005881
summary(TS[,2])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.001647 0.004798 0.005076 0.005053 0.005345 0.006746
summary(TS[,3])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.01677 0.04595 0.05292 0.05378 0.06079 0.11102
summary(TS[,4])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.09732 0.28162 0.38112 0.44060 0.53095 3.30433
summary(TS[,5])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.003112 0.003835 0.004003 0.004012 0.004178 0.005092
summary(TS[,6])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.003519 0.004254 0.004448 0.004460 0.004654 0.005750
summary(TS[,7])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.02251 0.03886 0.04480 0.04600 0.05173 0.10463
summary(TS[,8])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0568 0.2079 0.2988 0.3565 0.4359 2.9008
summary(TS[,9])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.384 1.918 2.075 2.092 2.242 3.502
summary(TS[,10])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.345 1.912 2.076 2.093 2.252 3.652
summary(TS[,11])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.5934 1.5712 2.0222 2.1859 2.6028 11.5787
summary(TS[,12])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1298 0.8929 1.5883 2.6582 2.9823 112.5042
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.995457 -1.994947 -1.946224 -1.559399 -1.995988 -1.99554 -1.954004
## [,8] [,9] [,10] [,11] [,12]
## [1,] -1.643466 0.09225532 0.09291883 0.1858899 0.6582028
EQMTS
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 3.981848 3.979812 3.787926 2.490755 3.983967 3.982179 3.818232 2.752803
## [,9] [,10] [,11] [,12]
## [1,] 0.0681228 0.07421034 0.8165807 14.36925
############################################################################## QUESTÃO 7
setwd("C:/Users/MarcelAlmeida/Cofco International/FSUSA - Coffee/Marcel/R/")
data_pop <- readxl::read_excel("pop2000.xlsx")
pop <- data_pop[,3]*1
#Características
K = 9
samplesz <- 50000+K*50
samplesz
## [1] 50450
subset <- subset(pop, pop >samplesz)
tt_subset <- matrix(0, ncol=3, nrow = 3) # estimates matrix
tt_subset[,1] <- mean(subset[,1])/(mean(subset[,1])-1) ### MM
tt_subset[,2] <- 1/mean(log(subset[,1])) ### MVM
tt_subset[,3] <- log(2)/log(median(subset[,1])) ### Mediana
#Calculando densidades
TT1_d <- density(tt_subset[,1]) #MM
TT2_d <- density(tt_subset[,2]) #MV
TT3_d <- density(tt_subset[,3]) #ME
#Plot
plot(TT1_d, xlim = c(-2,4), ylim = c(0, 5), lty = 1, lwd = 2, main = "Densidade dos Estimadores", xlab = " ", ylab = " ") # n1000
lines(TT2_d, lty = 2, lwd = 2, col = "red")
lines(TT3_d, lty = 3, lwd = 2, col = "blue")
legend("topright", legend = c("MM", "MV", "ME"), lty = 1:3, lwd = 2)
