1 Introducción

El Test de Jarque-Bera (JB) es una herramienta estadística esencial, desarrollada por Carlos Jarque y Anil K. Bera, que se utiliza para probar la hipótesis de normalidad de un conjunto de datos. Este test es particularmente significativo porque la normalidad es una suposición clave en muchos métodos estadísticos. La prueba evalúa si la distribución de un conjunto de datos difiere de la distribución normal en términos de asimetría (skewness) y curtosis(kurtosis).

El valor de la prueba compara la asimetría y la curtosis de la muestra con las expectativas de una distribución normal (donde la asimetría es cero y la curtosis es tres). Donde

  • \(H_0\): Los datos siguen una distribución normal
  • \(H_1\): Los datos no siguen una distribución normal

Los momentos muestrales son estadísticas clave que describen características esenciales de los datos, incluyendo la media, varianza, asimetría y curtosis. Estos momentos son cruciales en el Test de Jarque-Bera, el cual evalúa la normalidad de una distribución. La media indica la ubicación central de los datos y sirve como referencia para calcular la asimetría y curtosis. La varianza mide la dispersión de los datos alrededor de la media, esencial para entender su escala y dispersión. La asimetría (skewness) mide la simetría de la distribución; valores distintos de cero indican sesgo. Finalmente, la curtosis compara la “pesadez” de las colas de la distribución con una normal; valores superiores a 3 indican colas más pesadas (leptocúrtica), y menores que 3, colas más ligeras (platicúrtica). El Test de Jarque-Bera se basa en la asimetría y la curtosis para determinar desviaciones significativas de la normalidad, siendo un aspecto fundamental en diversas áreas de la estadística y la investigación.

1.1 Configuración inicial

La versión que se utilizará

version$version.string
## [1] "R version 4.3.2 (2023-10-31 ucrt)"

El preambulo contiene a la semilla y los paquetes utilizados. Al establecer una semilla set.seed(), se asegura que los resultados del código sean reproducibles y consistentes cada vez que se ejecutan.

#Preambulo
rm(list = ls())
options(scipen = 10) #notacion de los numeros grandes
set.seed(123)
library(ggplot2) #graficos
library(knitr) #cuadros

2 Simulación de Datos

En el análisis estadístico, la simulación de datos es una herramienta poderosa para comprender las propiedades y el comportamiento de diferentes distribuciones. En este contexto, vamos a generar dos conjuntos de datos simulados, cada uno con 100,000 observaciones, utilizando distribuciones uniforme y normal, respectivamente. Se elegírá una uniforme que coincida en media y varianza con al distribución estandar.Estas simulaciones nos permitirán explorar y visualizar las diferencias fundamentales entre estas dos distribuciones. Como ejemplo se usarán 2 distribuciones, normal estandár y la uniforme.

#numero de simulaciones
n <- 10^5

2.1 Muestra normal estándar

La distribución normal estándar es el caso particular de la distribución normal de media 0 (\(\mu\)=0) y varianza 1 (\(\sigma^2=1\)). La distribución normal es una de las más importantes en estadística debido a su presencia en muchos fenómenos naturales y teoremas estadísticos. Entonces \(X\sim N(0,1)\).

\[f_X(x)=\frac{1}{\sqrt{2\pi}}\exp\left\{ -\frac{x}{2} \right\}\]

#Distribuciones
mu <- 0
sigma2 <- 1
simulN <- rnorm(n, mu, sqrt(sigma2))

2.2 Muestra uniforme de media cero y varianza unitaria

Se usará la distribución uniforme

\[f_X(x)=\frac{1}{b-a}\] \[ E(X)=\frac{a+b}{2} \] \[V(X)=\frac{(b-a)^2}{12}\]

Se elegirá convenientemente los valroes de \(a\) y \(b\) para que cumplan \[\frac{a+b}{2}=\mu;\qquad \frac{(b-a)^{2}}{12}=\sigma^{2}\]

entonces \[ a = \mu - \sqrt{3\sigma^2} \]

\[ b = \mu + \sqrt{3\sigma^2} \]

#Distribucion uniforme
a <- mu - (3*sigma2)^0.5; a
## [1] -1.732051
b <- mu + (3*sigma2)^0.5; b
## [1] 1.732051
simulU <- runif(n, a, b)

2.3 Comparación de distribuciones

Se eligieron convenientemente los parámetros para que coincidan en media y varianza. También hemos generado una gráfica que compara las densidades muestrales.

#Tabla
estPob <- data.frame(
  media = c(mu,(a+b)/2),
  varianza = c(sigma2,(b-a)^2/12) 
)

et1 <- "N(0,1)";  et2 <- paste0("U(",round(a,2),",",round(b,2),")")

row.names(estPob) <- c(et1, et2)

kable(t(estPob))
N(0,1) U(-1.73,1.73)
media 0 0
varianza 1 1

Para complementar el análisis estadístico, hemos generado gráficos de densidad para ambas distribuciones. Estos gráficos ofrecen una representación visual clara de cómo se distribuyen los datos simulados en cada caso.

#Grafico
datos <- data.frame(1:n,simulU,simulN)

ggplot(data=datos) +
  geom_density(aes(x=simulU, color=et2))+
  geom_density(aes(x=simulN, color=et1)) +
  xlab("") + ylab("f(x)") + ggtitle("Densidades de las simulaciones") +
  theme_bw() +
  theme(legend.title = element_blank())

3 Momentos Muestrales Ordinarios y Centrales

3.1 Momentos muestrales ordinarios

Los momentos muestrales ordinarios ofrecen una visión significativa de la distribución de un conjunto de datos. Para una muestra dada de tamaño \(n\), el \(k\)-ésimo momento muestral ordinario se determina a través de la siguiente fórmula:

\[\hat{\mu}_{k}^\prime=\frac{1}{n} \sum_{i=1}^{n} x_{i}^{k} \]

Particularmente, el momento \(\hat{\mu}_{1}^\prime\) corresponde a la media muestral \(\bar{X}\).

La implementación en R para calcular estos momentos se presenta a continuación:

#Momentos muestrales ordinarios
momOrd <- function(x,k=1){
  n <- length(x)
  suma <- 0
  for (ii in 1:n) {
     suma <- suma+(x[ii])^k
  }
  return(suma/n)
}

Esta función, momOrd(), acepta un vector \(x\) y un entero \(k\) (con valor predeterminado 1), y retorna el \(k\)-ésimo momento muestral ordinario de x.

3.2 Momentos muestrales centrales

Los momentos muestrales centrales son cruciales para entender la variabilidad y la forma de una distribución alrededor de su media. Para una muestra de tamaño \(n\), el \(k\)-ésimo momento muestral central se calcula de la siguiente manera: \[\hat{\mu}_{k}=\frac{1}{n}\sum_{k=1}^{n}\left(x_{i}-\bar{X}\right)^{k}\] Es importante destacar que \(\hat{\mu}_{1}\) es igual a cero y \(\hat{\mu}_{2}\) es a la varianza sesgada, también como \(\hat{\sigma}^2\) (o asintóticamente insesgada).

La función de R para estos cálculos se muestra a continuación:

#Momentos muestrales centrales
momCent <- function(x,k=2){
  n <- length(x)
  media <- momOrd(x)
  suma <- 0
  for (ii in 1:n) {
    suma <- suma+(x[ii]-media)^k
  }
  return(suma/n)
}

Esta función, momCent(), toma un vector \(x\) y un entero \(k\) (por defecto 2, para calcular la varianza sesgada), y retorna el \(k\)-ésimo momento muestral central de \(x\). La función utiliza momOrd() para calcular la media de \(x\), que luego se usa en el cálculo del momento central.

3.3 Momentos muestrales de las simulaciones

En el cuadro siguiente, se ofrece un resumen comparativo entre los momentos muestrales y los momentos poblacionales. Dado que la muestra es considerablemente grande, se observa que los momentos muestrales se aproximan estrechamente a los poblacionales.

estMuest <- data.frame(
  c(n,et1,mu,momOrd(simulN)|>round(4),sigma2,momCent(simulN) |> round(4)),
  c(n,et2,(a+b)/2,momOrd(simulU)|>round(4),(b-a)^2/12,momCent(simulU)|> round(4))
)

colnames(estMuest) <- c("Normal","Uniforme")
row.names(estMuest) <- c("muestra","distribución","media poblacional",
                         "media muestral","varianza poblacional","varianza muestral")

kable(estMuest)
Normal Uniforme
muestra 100000 100000
distribución N(0,1) U(-1.73,1.73)
media poblacional 0 0
media muestral 0.001 0.0017
varianza poblacional 1 1
varianza muestral 0.9995 1.0066

4 Test de Jarque-Bera

El test de Jarque-Bera examina las propiedades de asimetría y curtosis de los datos para determinar cuán bien se ajustan a una distribución normal. La asimetría mide la desviación de los datos de la simetría perfecta, mientras que la curtosis evalúa la agudeza o planicidad de la distribución comparada con la curva normal.

El estadístico \(JB\) de la prueba es un indicador de la desviación conjunta de la asimetría y la curtosis respecto a aquellos valores esperados bajo una distribución normal (asimetría de cero y curtosis de tres). El estadístico se calcula de la siguiente forma:

\[JB=\frac{n}{6}\left(S^{2}+\frac{1}{4}(K-3)^{2}\right)\]

donde

  • \(n\) es el número de observaciones,

  • \(S\) es una medida de asimetría expresada como \[S=\frac{\hat{\mu}_{3}}{\hat{\mu}_{2}^{3/2}}=\frac{\frac{1}{n}\sum_{i=1}^{n}\left(x_{i}-\bar{X}\right)^{3}}{\left(\frac{1}{n}\sum_{i=1}^{n}\left(x_{i}-\bar{X}\right)^{2}\right)^{3/2}}\]

  • \(K\) es una medida de curtosis expresada como \[K=\frac{\hat{\mu}_{4}}{\hat{\mu}_{2}^{2}}=\frac{\frac{1}{n}\sum_{i=1}^{n}\left(x_{i}-\bar{X}\right)^{4}}{\left(\frac{1}{n}\sum_{i=1}^{n}\left(x_{i}-\bar{X}\right)^{2}\right)^{2}}\]

Jarque y Bera (1980) establecieron que, bajo la hipótesis nula de normalidad, el estadístico \(JB\) sigue asintóticamente una distribución chi-cuadrado con dos grados de libertad. El p-valor asociado con el estadístico JB se utiliza para decidir si rechazar o no la hipótesis nula de normalidad. Un p-valor pequeño (por lo general, menor a 0.05) indica que es improbable que la muestra provenga de una distribución normal.

El siguiente código en R implementa la función JarqueBera para calcular el estadístico JB y su p-valor asociado, permitiendo así realizar el test de Jarque-Bera en un conjunto de datos dado:

#Estadistico JB y su pvalue
JarqueBera <-function(x){
  n <- length(x)
  
  mu2 <- momCent(x,2)
  mu3 <- momCent(x,3)
  mu4 <- momCent(x,4)
  
  S <- mu3/(mu2^(3/2))
  K <- mu4/(mu2^2)
  
  JB <- (n/6)*(S^2+(K-3)^2/4)
  pvalue<-1-pchisq(JB,2)
  return(list("statJB"=JB,"pvalue"=pvalue, "skewness"=S, "kurtosis"=K-3))
}

Con esta función, puedes aplicar el test de Jarque-Bera a cualquier muestra de datos para evaluar la hipótesis de normalidad. Es importante realizar este test antes de proceder con métodos estadísticos que asumen normalidad, ya que una desviación significativa de esta suposición puede afectar la validez de los resultados y conclusiones del análisis.

4.1 Estadistico JB de las simulaciones

Una vez implementada la función JarqueBera en R, procedemos a aplicarla a nuestras simulaciones de datos. Recordemos que hemos generado dos conjuntos de datos: uno siguiendo una distribución normal \(N(0,1)\) y otro siguiendo una distribución uniforme \(U(-1.73, 1.73)\). El objetivo de aplicar el Test de Jarque-Bera a estas simulaciones es evaluar si los datos generados se ajustan a la distribución normal esperada, según los valores del estadístico JB y sus p-values asociados.

Los resultados obtenidos se muestran en la tabla siguiente:

#Jarque-Bera de las simulaciones
estJBN <- JarqueBera(simulN)
estJBU <- JarqueBera(simulU)
  
statJB <- data.frame(
  c(estJBN$skewness,estJBN$kurtosis,estJBN$statJB, estJBN$pvalue),
  c(estJBU$skewness,estJBU$kurtosis,estJBU$statJB, estJBU$pvalue)
)

colnames(statJB) <- c("Normal","Uniforme")
rownames(statJB) <- c("Asimetria","Exceso Curtosis","Estadistico JB","p-value")


kable(round(statJB,4))
Normal Uniforme
Asimetria -0.0039 0.0031
Exceso Curtosis -0.0099 -1.2049
Estadistico JB 0.6646 6049.0386
p-value 0.7173 0.0000

La interpretación de estos resultados es fundamental para entender la naturaleza de las distribuciones analizadas:

  • Distribución Normal (N(0,1)): El estadístico JB es relativamente bajo y el p-value es mayor que 0.05, lo que indica que no hay evidencia suficiente para rechazar la hipótesis nula de normalidad. Esto es coherente con nuestras expectativas, ya que los datos fueron generados a partir de una distribución normal.

  • Distribución Uniforme (U(-1.73, 1.73)): Por otro lado, observamos un valor muy alto para el estadístico JB y un p-value de 0.000 en la simulación de datos uniformes. Esto sugiere que hay evidencia significativa para rechazar la hipótesis nula de normalidad en este caso, lo cual era de esperar, dado que los datos provienen de una distribución uniforme, aunque provengan de una distribución simetrica,tiene propiedad de curtosis muy diferentes a las de una distribución normal.

5 Elipse de confianza del JB

A partir del valor crítico \(q\), que se elige como umbral para rechazar la hipótesis nula de normalidad. La condición que define la elipse es la siguiente

\[\frac{n}{6}\left(S^{2}+\frac{1}{4}(K-3)^{2}\right) = q\]

ordenando \[\frac{S^{2}}{6q\,/n}+\frac{(K-3)^{2}}{24q\,/n} = 1\] Donde \(S\) representa una medida de asimetría y \(K\) una medida de curtosis. Esta elipse está centrada en el origen, específicamente en \(S=0\) y \(K=3\), lo que corresponde a una distribución normal.La ecuación de la elipse centrada en el origen es \[\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}=1\]

Es decir, los semiejes de la elipse se calculan como \(a=\sqrt{6q/n}\) y \(b=\sqrt{24q/n}\)

La gráfica se puede automatizar

# Elipse de confianza
elipseConfianza <-function(conf,n){
  #conf: 1-alfa
  #n: numero de datos de la muestra
  
  q <- qchisq(conf, df = 2) #q critico o tabulado
  
  a <- sqrt(6 * q / n) 
  b <- sqrt(24 * q / n)
  
  Zona=paste0("Elipse de conf. (1-a =", conf, ", n =" ,format(n,scientific =T),")")
  
  # Ecuación paramétrica de la elipse 
  t <- seq(0, 2*pi, length.out = 1000)
  x <- a * cos(t) #Asimetria
  y <- b * sin(t) #Curtosis
  df <- data.frame(x, y)
  return(list(
    geom_path(data=df,aes(x,y,color=Zona),linewidth=1),
    labs(title = "Elipse de confianza",
         x = "Asimetría (S)",
         y = "Exceso de Curtosis (K-3)")
    ))
}

graficando

# Graficando la elipsis
conf <- 0.95 
q <- qchisq(conf, df = 2) #q critico o tabulado
n <- n #cantidad de la muestra
b <-  sqrt(24 * q / n)

# Limites del gráfico
lim_x <- 1.75 * b
lim_y <- 1.5 * b

#grafico
ggplot() + 
  elipseConfianza(conf,n) +   
  coord_fixed() +
  xlim(-lim_x, lim_x) +
  ylim(-lim_y, lim_y) +
  theme_bw()

5.1 Zona de rechazo y aceptación

Todos los puntos dentro de esta elipse se consideran distribuciones normales según el test de Jarque-Bera.

\[\frac{S^{2}}{6q\,/n}+\frac{(K-3)^{2}}{24q\,/n} \leq 1\]

y los que caen fuera rechazan la hipotesis nula de que los datos provienen de una normal \[\frac{S^{2}}{6q\,/n}+\frac{(K-3)^{2}}{24q\,/n} > 1\]

en el gráfico

# Elipse de confianza
conf <- 0.95
q <- qchisq(conf, df = 2) #q critico o tabulado
n <- n

a <- sqrt(6 * q / n) 
b <- sqrt(24 * q / n)

# Limites del gráfico
lim_x <- 1.75 * b
lim_y <- 1.5 * b

# Zona de rechazo y aceptación
x <- runif(500,-lim_x,lim_x)
y <- runif(500,-lim_y,lim_y)
df2 <- data.frame("y"=y,"x"=x)
df2$Zona <- ifelse((df2$x^2 / a^2 + df2$y^2 / b^2) <= 1, 
                   "Zona de Aceptación", "Zona de Rechazo")

# Graficar la elipse
ggplot() +
  geom_point(data=df2, aes(x,y,color=Zona)) +
  elipseConfianza(conf,n) +  
  theme_bw() +
  coord_fixed() +
  xlim(-lim_x, lim_x) +
  ylim(-lim_y, lim_y)

El área de una elipse es \[S=\pi \cdot a \cdot b\] en esta aplicación \[ S= \pi \cdot \sqrt{\frac{144q^2}{n^2}}=\frac{12\pi\cdot q}{n}\]

es decir, mientras aumente el nivel de confianza, aumentará el área, es decir habrá más tolerancia. Y mientras más aumente el tamaño de la muestra habrá menos tolerancia.

También podemos destacar que \[ \frac{a}{b}=\frac{\sqrt{6q\,/n}}{\sqrt{24q\,/n}}=\frac{1}{2} \]

la razón de los semiejes siempre será de 1 a 2 independimente del nivel de confianza (\(q\)) y el tamaño de muestra (\(n\)). Es decir, el test de Jarque-Bera tolera más la curtosis que la asimetría.

5.2 Elipse de confianza aplicado a las simulaciones

Para analizar y visualizar los resultados de las simulaciones, se ha utilizado la elipse de confianza del estadístico JB. En la tabla a continuación, se presentan los valores de asimetría (skewness) y curtosis (kurtosis) de las distribuciones simuladas, así como el tipo de distribución correspondiente.

df_simul <- data.frame(
  skewness = c(estJBN$skewness, estJBU$skewness),
  kurtosis = c(estJBN$kurtosis, estJBU$kurtosis),
  dist = c(et1,et2)
)

kable(df_simul)
skewness kurtosis dist
-0.0039011 -0.0099317 N(0,1)
0.0031371 -1.2048775 U(-1.73,1.73)

Posteriormente, se han graficado estos estadísticos en los ejes Simetria-curtosis. El primer y segundo gráfico son los mimos, pero a diferente escala.

ggplot() +
  elipseConfianza(conf,n) +  
  geom_point(data=df_simul,aes(x=skewness,y=kurtosis,color=dist)) +
  theme_bw() +
  coord_fixed() +
  xlim(-lim_x, lim_x) +
  ylim(-lim_y, lim_y)

ggplot() +
  elipseConfianza(conf,n) +  
  geom_point(data=df_simul,aes(x=skewness,y=kurtosis,color=dist)) +
  theme_bw() +
    xlim(-0.5, 0.5) +
  coord_fixed()

Vale la pena destacar que el estadístico JB de la distribución uniforme no se muestra en el primer gráfico debido a estar muy alejado. Por esto, en el segundo gráfico, podemos observar cuán alejado está de la zona de aceptación, lo que proporciona información importante sobre su ajuste a la distribución normal.

6 Aplicación: Test de Jarque Bera en N simulaciones

Lo que se ha abordado hasta ahora se limitó a una sola simulación de 100,000 datos con una distribución (la normal y uniforme), lo cual fue útil para propósitos explicativos y de comparación. Ahora, estamos preparando una aplicación que permitirá realizar N simulaciones y calcular las medidas promedio mediante una simulación de Montecarlo. Además, se proporcionará un gráfico que representará las N simulaciones en función de la cantidad de datos seleccionada. A continuación haremos un organigrama que nos ayudará a crear las funciones, una versión interactiva que requiere que el usuario ingrese datos. Y otra versión que no lo requiere.

6.1 Organigrama

En esta versión se controlo las fallas de las entradas como salidas del programa.

6.2 Explicación

En primer lugar, solicita al usuario que ingrese el número de simulaciones (N) y el tamaño de la muestra para cada simulación (n), asegurándose de que ambos valores sean válidos. Además, ofrece la opción de elegir entre dos distribuciones, “Normal” o “Uniforme”, validando la entrada del usuario.

Si se selecciona la distribución “Normal”, el programa solicitará al usuario que ingrese los parámetros de esta distribución, es decir, la media (mu) y la varianza (sigma2), validando que sean números adecuados. En el caso de la distribución “Uniforme”, se pedirá al usuario que ingrese los parámetros a y b, asegurando que a sea menor que b y que ambos sean valores numéricos válidos.

Además, el código permite al usuario decidir si desea generar gráficos para visualizar los resultados. En caso afirmativo, utiliza la biblioteca “ggplot2” para crear un gráfico que muestra las estadísticas de simetría y curtosis en una elipse de confianza, resaltando las áreas de aceptación y rechazo.

En el proceso de simulación, se inicializan vectores vacíos para almacenar los resultados de las estadísticas de Jarque-Bera (JB, pvalue, simetria, curtosis). Luego, se realizan las simulaciones correspondientes. Si se eligió la distribución “Normal”, se generan N muestras de la distribución normal con los parámetros especificados y se calculan las estadísticas de Jarque-Bera para cada muestra, almacenando los resultados en los vectores correspondientes. En el caso de la distribución “Uniforme”, se generan muestras de la distribución uniforme y se realiza el mismo proceso.

Una vez completadas las simulaciones, se crea un DataFrame (df) que contiene los resultados. Se calcula el promedio de las estadísticas de Jarque-Bera para todas las simulaciones y se muestra al usuario. Además, se determina el número de simulaciones donde se acepta la hipótesis nula (p-valor > 0.05) y donde se acepta la hipótesis alterna (p-valor <= 0.05), junto con las proporciones correspondientes.

6.3 Versión interactiva

#Version Interactiva
Nsimulaciones<-function(){
  # Solicitar el número de muestra de cada simulación, n
  repeat{
    cat("Introducir el número de simulaciones, N=")
    N <- readLines(n=1) 
    N <- as.numeric(N)
    if(N>=2){break}
    cat("Ingrese un número mayor o igual que 2","\n")
  }
  
  # Solicitar el número de muestra de cada simulación, n  
  repeat{
    cat("Introducir el número de muestra de cada simulación, n=")
    n <- readLines(n=1) 
    n <- as.numeric(n)
    if(n>=10){break}
    cat("Ingrese un número mayor o igual que 10","\n")
  }
  
  # Solicitar la elección de la distribución (Normal o Uniforme)
  repeat{
    cat("Elija una distribución: \n -Normal \n -Uniforme")
    dist <- readLines(n=1) 
    dist <- tolower(dist) #se acepta mayusculas o minusculas
    if(dist=="uniforme" | dist=="normal" ){break}
    cat("Ingrese un número mayor o igual que 10","\n")
  }
  
  # Solicitar los parámetros de la distribución Normal (si se elige)  
  if (dist=="normal"){
    cat("Elija los parametros de la normal: N(mu,sigma2)","\n")
    repeat{
      cat("mu=")
      mu <- readLines(n=1) 
      mu <- as.numeric(mu) 
      if(is.numeric(mu)){break}
      cat("Ingrese un número","\n")
    }

    repeat{
      cat("sigma2=")
      sigma2 <- readLines(n=1) 
      sigma2 <- as.numeric(sigma2) 
      if(sigma2>0){break}
      cat("Ingrese un número positivo","\n")
    }
  }
    
    
  # Solicitar los parámetros de la distribución Uniforme (si se elige)
  if (dist=="uniforme") {
    cat("Elija los parametros de la uniforme: U(a,b)","\n")
    repeat{
      cat("a=")
      a <- readLines(n=1) 
      a <- as.numeric(a) 
      if(is.numeric(a)){break}
      cat("Ingrese un número","\n")
    }
    
    repeat{
      cat("b=")
      b <- readLines(n=1) 
      b <- as.numeric(b) 
      if(b>a){break}
      cat("Ingrese un número mayor a ",a,"\n")
    }
    
  }
  
  # Solicitar si se desean ver gráficos (Sí o No)
  repeat{
    cat("Desea ver gráficos ¿Si o No?")
    graf <- readLines(n=1) 
    
    graf <- tolower(graf) #se acepta mayusculas o minusculas
    graf <- iconv(graf, to = "ASCII//TRANSLIT") #no acentos
    
    if(graf=="si" | graf=="no" ){break}
    cat("Ingrese un número mayor o igual que 10","\n")
  }
  
  # Inicializar vectores para almacenar resultados
  JB <- vector()
  pvalue <- vector()
  simetria <- vector()
  curtosis <- vector()
  
  # Realizar simulaciones
  if(dist=="normal"){
    for(ii in 1:N){
      muestra <-  rnorm(n, mu, sqrt(sigma2))
      est <- JarqueBera(muestra)
      JB <- c(JB,est$statJB)
      pvalue <- c(pvalue,est$pvalue)
      simetria <- c(simetria,est$skewness)
      curtosis <- c(curtosis,est$kurtosis)
    }
  }else{
    for(ii in 1:N){
      muestra <-  runif(n, a, b)
      est <- JarqueBera(muestra)
      JB <- c(JB,est$statJB)
      pvalue <- c(pvalue,est$pvalue)
      simetria <- c(simetria,est$skewness)
      curtosis <- c(curtosis,est$kurtosis)
    }    
  }
  
  # Crear un DataFrame con los resultados
  df <- data.frame(JB,pvalue,simetria,curtosis)
  
  
  # Mostrar el promedio de las estadísticas
  cat("Promedio de las", N, " simulaciones:", "\n")
  print(colMeans(df))
  
  # Contar cuántas simulaciones aceptaron o rechazaron la hipótesis nula
  N1 <- sum(df$pvalue > 0.05, na.rm = TRUE)
  N2 <- N-N1 
  
  cat("\n")
  cat("Simulaciones donde se acepto la hipotesis nula: ", N1 ,"\n")
  cat("Simulaciones donde se acepto la hipotesis alterna: ", N2 ,"\n")
  
  cat("Proporción donde se acepto la hipotesis nula: ", N1/N,"\n")
  cat("Proporción donde se acepto la hipotesis alterna: ", N2/N ,"\n")
  
  
  if(graf=="si"){
    if(!require("ggplot2", character.only = TRUE)){
      stop("Por favor installar el paquete ggplot2")}
    
    
    # Calcular valores para la elipse de confianza
    conf <- 0.95
    q <- qchisq(conf, df = 2) #q critico o tabulado
    n <- n
    
    a <- sqrt(6 * q / n) 
    b <- sqrt(24 * q / n)
    
    # Limites del gráfico
    lim_x <- 1.75 * b
    lim_y <- 1.5 * b
    
    x <- df$simetria
    y <- df$curtosis
    df2 <- data.frame("y"=y,"x"=x)
    df2$Zona <- ifelse((df2$x^2 / a^2 + df2$y^2 / b^2) <= 1, 
                       "Zona de Aceptación", "Zona de Rechazo")
    
    # Contar simulaciones fuera del gráfico
    fueraGrafico <- sum(x < -lim_x | x > lim_x |  y < -lim_y | y > lim_y)
    
    cat("\n")
    cat("Se realizo el gráfico con ", N, "Simulaciones","\n")
    cat("Nota: No se muestran ", fueraGrafico, 
        "simulaciones, por alejarse mucho de la zona de aceptación. ")
    
    ggplot() +
      geom_point(data=df2, aes(x,y,color=Zona)) +
      elipseConfianza(conf,n) +  
      theme_bw() +
      coord_fixed() +
      xlim(-lim_x, lim_x) +
      ylim(-lim_y, lim_y)
  }
}

6.4 Versión Función

Aqui no se pide ninguna entrada de usuario

#Version Funcion
NsimulacionesFuncion <-function(N,n,dist,parametro1, parametro2,graf){
  dist <- tolower(dist)
  graf <- tolower(graf) #se acepta mayusculas o minusculas
  graf <- iconv(graf, to = "ASCII//TRANSLIT") #no acentos
  
  if(N<2){stop("Ingresar N >= 2")}
  if(n<10){stop("Ingresar n >= 10")}
  if(!(dist=="normal"| dist=="uniforme")){stop("Distribuciones aceptadas: Normal o Uniforme")}
  if(dist=="normal" & parametro2 < 0){stop("La varianza es positiva")}
  if(dist=="uniforme" & parametro2 <= parametro1){stop("En una uniforme: b>a")}
  if(!(graf=="si"| graf=="no")){stop("Grafico: Si o No")}
  
  # Inicializar vectores para almacenar resultados
  JB <- vector()
  pvalue <- vector()
  simetria <- vector()
  curtosis <- vector()
  
  # Realizar simulaciones
  if(dist=="normal"){
    for(ii in 1:N){
      muestra <-  rnorm(n, parametro1, sqrt(parametro2))
      est <- JarqueBera(muestra)
      JB <- c(JB,est$statJB)
      pvalue <- c(pvalue,est$pvalue)
      simetria <- c(simetria,est$skewness)
      curtosis <- c(curtosis,est$kurtosis)
    }
  }else{
    for(ii in 1:N){
      muestra <-  runif(n, parametro1, parametro2)
      est <- JarqueBera(muestra)
      JB <- c(JB,est$statJB)
      pvalue <- c(pvalue,est$pvalue)
      simetria <- c(simetria,est$skewness)
      curtosis <- c(curtosis,est$kurtosis)
    }    
  }
  
  # Crear un DataFrame con los resultados
  df <- data.frame(JB,pvalue,simetria,curtosis)
  
  
  # Mostrar el promedio de las estadísticas
  cat("Promedio de las", N, " simulaciones:", "\n")
  print(colMeans(df))
  
  # Contar cuántas simulaciones aceptaron o rechazaron la hipótesis nula
  N1 <- sum(df$pvalue > 0.05, na.rm = TRUE)
  N2 <- N-N1 
  
  cat("\n")
  cat("Simulaciones donde se acepto la hipotesis nula: ", N1 ,"\n")
  cat("Simulaciones donde se acepto la hipotesis alterna: ", N2 ,"\n")
  
  cat("Proporción donde se acepto la hipotesis nula: ", N1/N,"\n")
  cat("Proporción donde se acepto la hipotesis alterna: ", N2/N ,"\n")
  
  
  if(graf=="si"){
    if(!require("ggplot2", character.only = TRUE)){stop("Por favor installar el paquete ggplot2")}
    
    
    # Calcular valores para la elipse de confianza
    conf <- 0.95
    q <- qchisq(conf, df = 2) #q critico o tabulado
    n <- n
    
    a <- sqrt(6 * q / n) 
    b <- sqrt(24 * q / n)
    
    # Limites del gráfico
    lim_x <- 1.75 * b
    lim_y <- 1.5 * b
    
    x <- df$simetria
    y <- df$curtosis
    df2 <- data.frame("y"=y,"x"=x)
    df2$Zona <- ifelse((df2$x^2 / a^2 + df2$y^2 / b^2) <= 1, 
                       "Zona de Aceptación", "Zona de Rechazo")
    
    # Contar simulaciones fuera del gráfico
    fueraGrafico <- sum(x < -lim_x | x > lim_x |  y < -lim_y | y > lim_y)
    
    cat("\n")
    cat("Se realizo el gráfico con ", N, "Simulaciones","\n")
    cat("Nota: No se muestran ", fueraGrafico, "simulaciones, por alejarse mucho de la zona de aceptación. ")
    
    ggplot() +
      geom_point(data=df2, aes(x,y,color=Zona)) +
      elipseConfianza(conf,n) +  
      theme_bw() +
      coord_fixed() +
      xlim(-lim_x, lim_x) +
      ylim(-lim_y, lim_y)
  }
}

6.5 Ejemplos de la versión no interactiva

#Ejemplos
NsimulacionesFuncion(100,100,"uniforme",5,10,"Si")
## Promedio de las 100  simulaciones: 
##          JB      pvalue    simetria    curtosis 
##  6.10429291  0.05535880 -0.01609664 -1.17386152 
## 
## Simulaciones donde se acepto la hipotesis nula:  42 
## Simulaciones donde se acepto la hipotesis alterna:  58 
## Proporción donde se acepto la hipotesis nula:  0.42 
## Proporción donde se acepto la hipotesis alterna:  0.58 
## 
## Se realizo el gráfico con  100 Simulaciones 
## Nota: No se muestran  0 simulaciones, por alejarse mucho de la zona de aceptación.

NsimulacionesFuncion(1000,1000,"normal",0,1,"Si")
## Promedio de las 1000  simulaciones: 
##             JB         pvalue       simetria       curtosis 
##  2.07043732542  0.49151049034 -0.00003058555 -0.00218536352 
## 
## Simulaciones donde se acepto la hipotesis nula:  952 
## Simulaciones donde se acepto la hipotesis alterna:  48 
## Proporción donde se acepto la hipotesis nula:  0.952 
## Proporción donde se acepto la hipotesis alterna:  0.048 
## 
## Se realizo el gráfico con  1000 Simulaciones 
## Nota: No se muestran  2 simulaciones, por alejarse mucho de la zona de aceptación.
## Warning: Removed 2 rows containing missing values (`geom_point()`).

NsimulacionesFuncion(1000,1000,"normal",-5,5,"Si")
## Promedio de las 1000  simulaciones: 
##           JB       pvalue     simetria     curtosis 
##  2.018014319  0.506122724 -0.001063820 -0.004666395 
## 
## Simulaciones donde se acepto la hipotesis nula:  939 
## Simulaciones donde se acepto la hipotesis alterna:  61 
## Proporción donde se acepto la hipotesis nula:  0.939 
## Proporción donde se acepto la hipotesis alterna:  0.061 
## 
## Se realizo el gráfico con  1000 Simulaciones 
## Nota: No se muestran  2 simulaciones, por alejarse mucho de la zona de aceptación.
## Warning: Removed 2 rows containing missing values (`geom_point()`).

NsimulacionesFuncion(1000,1000,"uniforme",0.5,3,"No")
## Promedio de las 1000  simulaciones: 
##                     JB                 pvalue               simetria 
## 60.3324215788171116515  0.0000000000003882831  0.0017739433448172239 
##               curtosis 
## -1.1992989600911279080 
## 
## Simulaciones donde se acepto la hipotesis nula:  0 
## Simulaciones donde se acepto la hipotesis alterna:  1000 
## Proporción donde se acepto la hipotesis nula:  0 
## Proporción donde se acepto la hipotesis alterna:  1
NsimulacionesFuncion(100,15000,"normal",100,500,"No")
## Promedio de las 100  simulaciones: 
##             JB         pvalue       simetria       curtosis 
##  2.05853773623  0.48806317163  0.00001998704 -0.00086187137 
## 
## Simulaciones donde se acepto la hipotesis nula:  94 
## Simulaciones donde se acepto la hipotesis alterna:  6 
## Proporción donde se acepto la hipotesis nula:  0.94 
## Proporción donde se acepto la hipotesis alterna:  0.06

  1. Economista. Estudiante del Master Ciencias Actuariales de la Universidad de Barcelona, correo: ↩︎