Este borrador .html contiene la segunda parte del curso de introducción a R de los programas de Ph.D. y Master en Economia y Finanza Computacional de UHU-UNIA. Este curso es autocontenido, de todas formas se aconseja la reproducción paso por paso del codigo de este fichero, siendo presentes interrupciones de tipo didactico condicionadas a la explicación de los conceptos del curso.
Podriamos no repetirlo, pero aquí abajo hay dos ejemplos (el primero en entorno Mac, el segundo en Widnows) de declaración de la carpeta de trabajo. Esto paso se hace encesario cuando se queire trabajar con datos que requieren el utilizo del disco duro del ordenador
setwd("/Users/nicola/Desktop/train")
setwd("C:/Users/nicola/OneDrive/Desktop/train")Creamos dos vectores sorteando valores de una distribución normal estandardizada:
x <- rnorm(50)
print(x)
y <- rnorm(x)
print(x)en el caso de x, hemos pedido 50 elementos de una distribucion \(N(0,1)\). para. crear y, hemos pedido un vector identicamente normalizado cogiendo aleatoriamente elementos que esten presentes en x.
Vamos a mirarnos un sencillo scatterplot:
plot(x, y)hechamos un vistazo al global environment:
ls()o, el que da igual y es ecuivalente a la función “ls()”:
objects()ahora, como hicimos arriba para quitar los objetos de la sesión global, utilizamos la función remove:
rm(x,y)creamos una vecror que va de 1 a 20:
x <- 1:20hay que notar que esto es distinto del comando concatenate:
x <- c(1,20)
x <- 1:20Vamos creando un vector de desviaciones en función de x:
z<-1:20
w<- 1+sqrt(x)/2Como ya sabemos, el dataframe es el conjunto de objetos utilizados en la analysis. vamos declarando sus contendidos:
ourdata <-data.frame(x=x, y=x+rnorm(x)*w)
dummy <-data.frame(x=x, y=x+rnorm(x)*w)ahora x será igual a x (covariata), y será igual a una series de sucesos de x multiplicados por su desviacion estandar. (para trabajar mas comodamente, a veces vale la pena crear una copia inmediata del data frame, un data frame “dummy”)
print(dummy)############################################
para la tilde, option+ñ en teclado mac. para x.186,
Se trata, sencillamente, de llamar la función “lm”: [^1]
[^]: Para la tilde, option+ñ y luego la tecla espacio en teclado mac. Para x.186, ctrl+ñ y luego el tecaldo espacio.
lineal1 <- lm (y~x, data=ourdata)
summary(lineal1)y dado que conocemos las desviaciones estandar, podemos hacer una regresion pesada:
linealw <- lm(y~x, data=ourdata, weight =1/w^2)
summary(linealw)el mando attach sirve para que los vectores columnas en el dataframe sean identificables como variables esto nos permite evitar de declarar la opción lm(data=)
attach(dummy)
view(dummy)lineal1 <- lm (y~x)
summary(lineal1)las suites basicas de R nos ponen a disposición otra herramientas utiles. SI quisieramos hacer un regresión lineal no parametrica por ejemplo:
frl <-lowess(x,y)
summary(frl)plot(x,y)vamos visualizando la linea de regresión no parametrica:
plot(x, frl$y)ahora visualizamos la linea de las y:
plot(x, lineal1$y)ahora la verdadera linia de regresión ipotetica con intecepta 0 e pendiente 1:
abline(0,1, lty=3)La linea de regresión lineal:
abline(coef(lineal1))la linea de regresión pesada:
abline(coef(linealw))y para diferenciar la regresión pesada de la normal, utilizar el color verde:
abline(coef(linealw), col ="green")Podemos quitar el dataframe de la. directory de busqueda:
detach()así que:
linealwrong <- lm(y~x)nos dará un mensaje de error
Vamos dibujando un grafico de analisis mas elegante:
plot(fitted(lineal1), resid(lineal1), xlab="fitted values", ylab= "Residuals", main= "Residuals vs Fitted")el grafico de arriba nos resultaría util para verificar si hay heteroskedasticidad en los residuos.
Verificamos ahora que el tercero y cuarto momento de la distribución sean normales:
qqnorm(resid(lineal1), main ="residuals RankitPlot")para acabar, vamos limpiando todo el entorno:
rm(dummy,lineal1, linealw, ourdata, x, frl, w)############################################################
Dejando a parte las aprlicaciones de corte transversal o de series temporales, muy a menudo las aplicaciones parametricas se basan en datos a dos dimensiones. Los datos de panel poseen un identificador (dimensión de corte transversal) y una tendencia (dimensión temporal). Antes de todo, instalamos la libreria para modelos lineales de panel:
library("plm")vais a notar que no hay nada que se llame así en libreria. Por lo tanto instalamos el paquete del CRAN:
install.packages("plm")
library(plm)Sacamos ahora lo tanto los datos Grunfeld, un panel de empresas estadounidenses de los años 30 muy utilizados para ejemplos econometricos introductorios de panel. La frecuencia de los datos es anual, y las variables de interes son las inversiones (la formación de capital fijo bruto), el valor de mercado (la cantitad nominal de capitalización) y el capital fisico a disposición (los assets de cada empresa). El identificador es “firm”, la tendencia es “year”.
data("Grunfeld")hay que subrayar que el dataset Grunfeld, que carga un data frame en la memoria volatíl del Global Environment de R pertenece al paquete plm. Para sacar los datos sin necesariamente activar plm en libreria, se podría hacer:
data("Grunfeld", package="plm")head(Grunfeld)
View(Grunfeld)
Print(Grunfeld)notar arriba que print no ha funcionado. porqué?1
Proximo step: Declarar los indices que identifican el panel de datos. Es decir, decalrar las dos dimensiones:
Grunfeldpanel<-pdata.frame(Grunfeld, index=c("firm", "year"), drop.index=TRUE, row.names=TRUE)
head(Grunfeldpanel)Que ha pasado? pues ahora tenemos las dos dimensiones delp panel identificadas. Hemos dicho a R de quitar los vectores de las dimensiones, y de ubicar el indice a lado del conjunto de datos.
head(attr(Grunfeldpanel, "index"))el objeto panel data frame es mas complejo del data frame: conteine facotres debidos a la prensencia de los indices. Ademas, hay metodos especificos, como summary() o as.matrix(), que reaccionan de forma contextual cuando se aplican a un paquete como plm.
Summary() es una analysis de la varianza interdimensional: nos da una idea de la suma del cuadrado del error de cada variable, el porcentaje de variación debido a la dimension de corte trnasversal o a la dimensión temporal, y ademas nos restituye unas estadisticas descriptivas:
summary(Grunfeldpanel)si quisieramos analizar solo y exclusivamente el capital, cual sería la syntaxis correcta?
summary(Grunfeldpanel$capital)as.matrix(), al contario, nos expressa el panel en el formado de matriz equivalente:
head(as.matrix(Grunfeldpanel$capital))A veces, para crear modelos dinamicos, es necesario considerar una serie con un periodo de retraso. Eso Equivale a crear un vector equivalente a la variable original que esté depslazado hacia delante. Pongamos que quiera obtener un retardo de la variable capital:
head(lag(Grunfeldpanel$capital, 1))Como obtener dos retardos?
head(lag(Grunfeldpanel$capital, 0:2))A veces, necesitamos trabajar con las variaciones, en lugar de los niveles. O con las tasas de variacion. para hacer eso, hay que calcular una diferencia. una diferencia de primer orden por ejemplo se calcularía con:
head(diff(Grunfeldpanel$capital, 1))una del segundo orden con:
head(diff(Grunfeldpanel$capital, 2))para las dos:
head(diff(Grunfeldpanel$capital, 0:2))hay claramente otras posibilidades. Estando nosotros en un entorno de Panel de datos, es posible efectuar las transformaciones necesarias para calcular la variacion diferencial dentro cada corte transversal y entre los cortes (las transformaciones necesarias para los estimadores within y between):
head(Within(Grunfeldpanel$capital))
head(between(Grunfeldpanel$capital))Pasamos a la estimación. En el entorno de datos de panel, tenemos:
la formulación sigue la sintaxis que hemos vistgo para el modelo lineal. Efectivamente, un modelo pooled en panel de datos es equivalente a un modelo lineal de corte transversal (y, por supuesto, un modelo a efectos fijos es equivalente a un modelo de corte transversal con un indicador categorico)
G.fe <- plm(inv~value+capital, data = Grunfeld, model = "within")
G.po <-plm(inv~value+capital, data=Grunfeld, model="pooling")
G.re <- plm(inv~value+capital, data = Grunfeld, model = "random")claramente, siendo plm() un metodo con sus propias listas y sus propias funciones, los inputs que damos serán posicionales. Si de todas formas no tenemos clara la posición de cada objeto de la lista, podemos utilizar la etiqueta del objeto/opción de la función.2
G.po <-plm(inv~value+capital, Grunfeld, model="pooling")
G.po <-plm(inv~value+capital, Grunfeld, "pooling")Los valores resultantes de la estimación se guardarán en los objetos listas que acabamos de nombrar. Para visualizar el resultado de la estimacion sin utilizar la notación “$”, vamos a por un summary()
summary(G.fe)
summary(G.po)
summary(G.re)PARA UNA OVERVIEW SUPER-EXHAUSTIVA DE LOS MODELOS DE PANEL DE DATOS EN R: https://cran.r-project.org/web/packages/plm/vignettes/A_plmPackage.html
#######################################################################
Una serie historica es un vector de datos ordenados temporalmente. Creamos un vector de 10 valores:
Tsexample<-c(1:10)Para utilizar las herramientas de series historicas de R, hay que cambiar de objeto! Utilizando nuevamente la funzion as.(), vamos a declarar nuestra serie:
tsdeclared<-as.ts(Tsexample)intentamos hacerlo mejor!
tsdeclaredbetter<-as.ts(Tsexample, start=1991, frequency=1)
print(tsdeclaredbetter)Vamos a utilizar una base ya lista. Los datos que vamos a utilizar representa el flujo anual del rio Nilo durante 100 años, medido en 108m3:
data("Nile")
help("Nile")Funciones como print(), length(), head() o tail() proveen, como siempre, informaciones utiles sobre nuestra serie historica:
print(Nile)
length(Nile)
head(Nile, 10)
tail(Nile, 5)La sub-seleccion con las corchetas cuadradas tambien funciona: para imprimir el valor numero 2 de la serie historica, por ejemplo:
print(Nile[2])En este contexto, la función plot() reconoce el objeto como serie historica y prepara un index temporal en el eje de las x:
plot(Nile)como antes, podemos modificar unas cuantas opciones:
plot(Nile, xlab = "Año", ylab = "Vol, 100 de ML. de M3")como ya vimos, con la opcion main podemos nombrar el grafico. Ademas, podemos seleccionar distintos estilos:
plot(Nile, xlab = "A?o", ylab = "Vol, 100 de ML. de M3", mai="Volumen anual del flujo del Nilo en Abu-Simbel", type="b")Las funciones start() y end() nos dan el indice de serie de la primera y ultima observacion de la serie.
print(start(Nile))
print(end(Nile))las funciones deltat() restituye el numero de obsrvaciones por cada interval de tiempo la funcion frequency() el numero de observaciones por unidad de tiempo (claramente en este caso son parecidos!)
print(deltat(Nile))
print(frequency(Nile))Asì como vimos con el ejemplo de panel de datos, todas las herramientas de transformación de variable siguen siendo las mismas!
logNile<-log(Nile)
diffNile<-diff(Nile,1)
diff2Nile<-diff(Nile,2)plot(logNile)
plot(diffNile)El modelo de ruido blanco es un tipico modelo de series historicas. Se trata de un modelo estacionario, un caso peculiar de modelo autoregresivo integrado de media movil (ARIMA). El modelo posee un promedio y una desviación constantes en el tiempo, y ningun elemento autoregresivo o de media movil ARIMA(0,0,0). Intentamos simularlo:
ruidoblanco<- arima.sim(model = list(order = c(0, 0, 0)), n = 100)por defecto, hemos obtenido un modelo de ruido con promedio 0 y desviación estandar 1!
ts.plot(ruidoblanco)ruidoblanco2<-arima.sim(model = list(order=c(0, 0 ,0)), n=100, mean = 50, sd=5)
ts.plot(ruidoblanco2)Estamos convencido que el modelo aleatorio de arriba se adecue a los datos del Nilo. Que tal comprovar la bontad de ajuste con una estimación?
Arimanilo<-arima(Nile, order=c(0,0,0))
summary(Arimanilo)
print(Arimanilo)Ahora podemos calcular comodamente promedio, varianza y desviaci?n estandar de las serie:
print(mean(Nile))
print(var(Nile))
print(sd(Nile))Imaginamos ahora que el flujo del Nilo, sus valores contemporaneos, dependan de la historia pasada del rio, de forma no estacionaria. Un modelo historico no estacionario tipico es el paseo aleatorio:
paseoaleatorio <- arima.sim(model = list(order = c(0, 1, 0)), n = 50)
ts.plot(paseoaleatorio)la diferencia de primer orden de un paseo aleatorio es estacionaria por definicion:
ts.plot(diff(paseoaleatorio))una categoria interesante de paseo aleatorio es el paseo aleatorio con deriva fundamentalmente es un paseo aleatorio con un promedio distinto de 0:
paseoaleatorioderiva <- arima.sim(model = list(order = c(0, 1, 0)), n = 50, mean=5)
ts.plot(paseoaleatorioderiva)Arimanilorw<-arima(Nile, order=c(0,1,0))
summary(Arimanilorw)
print(Arimanilorw)claramente, examinando el AIC entre Arimanilo y Arimanilorw podemos escoger el modelo mas adecuado:
print(Arimanilo)
print(Arimanilorw)Analizamos los modelos simulados juntamente en un grafico:
plot.ts(cbind(paseoaleatorio, paseoaleatorioderiva, ruidoblanco))Para una analisis mas exaustiva de las herramientas de series historicas: https://rpubs.com/odenipinedo/time-series-analysis-in-R https://cran.r-project.org/web/views/TimeSeries.html
#######################################################################
En programación, queremos automatizar determinados procesos para que se repitan constantemente. Los tipos de bucles a nuestra desposición son normalmente de “for”, “while” y de “repeat”.
Imaginamos crear un objeto contenente un escalar 0:
x_por <- 0utilizaremos este objeto como base del bucle. Ahora necesitamos un objeto mas como elemento de iteración: suponemos querer añadir un escalar de valor 1, de forma progresiva, a nuestro vector 0. Primero, necesitamos la intestación de loop:
for(n in 1:10) {
x_por <- x_por+1
print(x_por)
}
rm(x_por)de forma equivalente, podemos crear nosotros mismos el elemento de iteración:
x_por <- 0
x_iter<-(1:10)
for(i in x_iter) {
x_por<-x_por+1
print(x_por)
}claramente, las dos formulaciones de arriba nos darán un resultato distinto de:
x_por <- 0
x_iter<-(1:10)
for(i in x_iter) {
x_por<-x_por+i
print(x_por)
}DOnde la progresión geometrica es de razón n(t+1)=n(t)+n(t-1)+1
Los bucles con la condición “mientras” expresan una condición temporanea de ejecución del contenido del bucle.
Como antes, vamos creando la base del bucle:
x_mientras <-0y ahora construimos el bucle condicional:
while(x_mientras<10) {
x_mientras<-x_mientras+1
print(x_mientras)
}Esto nos permete ejecutar la operación hasta que la suma del valor de la base del bucle con un escalar 1 no llegue al limite de la condición de tener que ser menor que 10.
Los Bucles con repeat siguen la misma logica de los que vimos antes, pero nos permiten introducir una condición de interrupción de la procedura en curso expressando una condición adicional con el operador if() dentro del cuerpo principal del bucle.
x_repetir<-0
repeat{
x_repetir<-x_repetir+1
print(x_repetir)
if(x_repetir>10) {
break
}
}Notar como el algoritmo recursivo haya considerdado los valores hasta 11. Para que se pare en 10, claramente hay que modificarr el argumento del operador if de la función de “break” de la operación con: if(x_repetir>+10).3
#######################################################################
EJEMPLO 1
Evimatos de definir la base del loop y la dejamos en el cuerpo principal del bucle:
for(i in 1:10) {
x1 <- i^3
print(x1)
}
De todas formas, aquí hay un problema. En lugar de print, querríamos guardar los valores e iteración en un vector. Por lo tanto, vamos creando un vector, le ponemos valor nulo, y aplicamos con la función concadenar la progresión cubica que nos interesa.
xcubo <- numeric()
for(g in 1:10) {
xcubo <- c(xcubo, g^3)
print(xcubo)
}
print(xcubo)Notar como el print del bucle considera toda la progresión, el print externo al bucle solo el resultado final!
EJEMPLO 2: BUCLE DE UN BUCLE CON FOR
Aprovechamos para ampliar aún mas nuestro conocimiento de los bucles montando un doble bucle de por. Creando dos indices i e j de 1 a 5 y de 1 a 3, nuestro bucle quiere que se peguen las primeras 5 minusculas del alfabeto con las primeras 3 MAYUSCULAS:
xletras<-character()
for(i in 1:5) {
for(j in 1:3) {
xletras <- c(xletras, paste(letters[i], LETTERS[j], sep = "_"))
}
}
print(xletras)de la misma forma, porqué no multiplicar los primeros cinco numeros reales \((n=[1:5]\in R)\) reales por los primeros 3 \((n=[1:3]\in R)\)?
xnum<-numeric()
enteros<-integer(5)
enteros<-c(1:5)
for(i in 1:5) {
for(j in 1:3) {
xnum <- c(xnum, paste(enteros[i]*enteros[j]))
}
}
print(xnum)
print(integer(1:5))EJEMPLO 3: INTERRUPCIÓN DE BUCLE
Pongamos de querer calcular el cuadrado de los primero 10 numeros enteros, menos el 3, el 6 y el 7. Con %in% vamos a subrayar la excepcion, con “next” indicamos a R que hay que no hay que considerar esos valores en el mando sucesivo.
for(i in 1:10) {
if(i %in% c(3, 6, 7)) {
next
}
xcuad <- i^2
print(xcuad)
}Como hacer coger, al contrario, los valores cuadrados que no sean los de 3, 6 y 7 sin modificar excesivamente el codigo? Con el simbolo de negación, “!”.
for(i in 1:10) {
if(i %!in% c(3, 6, 7)) {
next
}
xcuad <- i^2
print(xcuad)
}hay algo raro……. el simbolo no parece existir….habrá que crear esa función utilizando unas que ya están reconocidas por R!
`%!in%` = Negate(`%in%`)Gracias a esta nueva función, la modifica del bucle resulta sencillisima:
for(i in 1:10) {
if(i %!in% c(3, 6, 7)) {
next
}
xcuad <- i^2
print(xcuad)
}Para una overview muy completa sobre lso bucles: https://statisticsglobe.com/loops-in-r/
############################################################
data(trees)
help(trees)Volvemos a un ejemplo real para mostrar un bucle de función lineal. El dataset “trees” contiene diametro, altura y volumen de 31 arboles de cereza. Nuestra variable objetivo será por lo tanto el diametro de los arboles, y nuestra idea sería efectuar una serie de regresiones lineales añadiendo progresivamente todas las covariatas posibles (justo dos, en este caso: altura y volumen)
summary(trees)creamos un data.frame de preuba, treesnew:
treesnew<-treesCrear una lista que luego nos permita alamcenar los resultados de los summaries repetidos que queremos:
mod_summaries <- list() Preparar la intestación del loop con todos los valores de la segunda columna hasta la ultima (el conjunto de los regresores):
for(i in 2:ncol(treesnew)) { crear un vector con los nombres de los regressores de la segunda a la i-esima (tercera y ultima) columna:
predictors_i <- colnames(treesnew)[2:i] indicar el metodo de almacenamiento de los resultados de las regresiones:
mod_summaries[[i - 1]] <- summary(lm(Girth ~ ., treesnew[ , c("Girth", predictors_i)]))
}
print(mod_summaries)
Esto concluye el curso.
Porqué no intentais cambiar la P mayuscola de Print con una minuscola?↩︎
El output sigue siendo igual. Si quitamos formula y data sigue funcionando. Si le quitamos model no, porqué tras formula y data hay otro objeto y no se trata de model!!!↩︎
El mando break es utilizable en cualquier tipologia de bucle↩︎