Machine Learning to identify the best factor models for portfolio design

Problem Setup

En el siguiente trabajo se presenta la resolución, elaboración e implementación de los mejores algoritmos para diseñar carteras financieras para los clientes de SmartPort. En SmartPort, los métodos para seleccionar acciones y diseñar carteras han sido el análisis fundamental tradicional utilizando razones financieras simples y análisis técnico básico. No obstante, en búsqueda de incrementar la confiabilidad de sekección de activos y brindar mayor soporte estadístico a la prueba, el trabajo implementa algoritmos de aprendizaje automático para una selección más efectiva en el mercado de acciones estadounidenses.

Análisis de datos

Tipos de datos para el análisis financiero

Para llevar a cabo un adecuado análisis financiero es importante obtener una lectura de los datos históricos reportados por las empresas, industrias o mercados en cuestión, ya que mediante el análisis de estos podremos evaluar el desempeño histórico y bien, generar un panorama real de los posibles futuros resultados.

En Finanzas, existen los siguientes 3 principales tipos de datos:

  1. Series de tiempo (Time-Series): una serie temporal o cronológica es una sucesión de datos medidos en determinados momentos y ordenados cronológicamente. Los datos pueden estar espaciados a intervalos iguales o desiguales. En esta estructura puede haber muchos períodos e información para un tema (empresa, industria, mercado, índice, etc.)

  2. Estructura transversal (Cross-sectional structure): los datos transversales o una sección transversal, en estadística y econometría son un tipo de datos recopilados mediante la observación de muchos sujetos (como individuos, empresas, países o regiones) en un punto o período de tiempo. En esta estructura se suelen tener muchos “sujetos”, para un solo período.

  3. Estructura de datos de panel (Panel-data structure): los datos de panel estructurados, son datos que se derivan de un número (generalmente pequeño) de observaciones a lo largo del tiempo en un número (generalmente grande) de unidades transversales como individuos, hogares, empresas o gobiernos. Es una combinación de series de tiempo y datos transversales. En esta estructura, tenemos más de una asignatura, y para cada asignatura podemos tener más de un período.

Bases de datos utilizadas para el análisis

Para la resolución es esta situación problema, trabajaremos con un panel de datos de firmas públicas estadounidenses reales. Puntualmente, trabajaremos con las siguientes 3 bases de datos:

  1. usfirms.xlsx: lista de firmas públicas estadounidenses con información general de cada firma.
  2. us2020q.xlsx / data: panel de datos con datos financieros anuales históricos para todas las empresas públicas de EE. UU.
  3. us2020q.xlsx / data dictionary: diccionario para cada variable en el conjunto de datos del panel.
# Descargamos las bases de datos desde internet.

options(scipen=999)

library(readxl)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(PerformanceAnalytics)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
download.file("http://www.apradie.com/datos/us2020q.xlsx",
              "us2020q.xlsx", mode="wb")

download.file("http://www.apradie.com/datos/usfirms.xlsx",
              "usfirms.xlsx", mode="wb")

Importamos los datos trimestrales de panel de todas las empresas estadounidenses:

uspanel <- read_excel("us2020q.xlsx",sheet = "data")

Importamos el diccionario de datos para las variables de este conjunto de datos:

dictionary <- read_excel("us2020q.xlsx",sheet = "data dictionary")

Importamos la lista de todas las empresas públicas estadounidenses:

usfirms <- read_excel("usfirms.xlsx")

Pasos realizados de manejo de datos (data merging, data selection, cálculos de variables)

Creo un conjunto de datos con solo las columnas relevantes para fusionar, y a su vez, hago la combinación indicando que coincida con ambos archivos usando la columna firm que se encuentra en ambos conjuntos de datos.

usfirms1 <- usfirms %>% select(firm, CompanyName, status, naics1)
uspanel <- merge(uspanel, usfirms1, by="firm")    

uspanel

Para la elaboración de este trabajo se analizó única y exclusivamente las empresas pertenecientes a la industria de manufactura de los Estados Unidos, a continuación creamos dos diversos paneles:

uspanel_2020 <- uspanel %>%
          filter(status=="active", year==2020, fiscalmonth==12, naics1 == "Manufacturing")

uspanel_2020
#Este primer panel ha sido creado filtrando el uspanel por status, año, mes fiscal e industria; nos servirá para posteriormente hacer un breve análisis descriptivo de la industria al año 2020. 

uspanel <- uspanel %>%
          filter(status=="active", naics1 == "Manufacturing")

uspanel
#Este segundo panel ha sido creado para el desarrollo del trabajo en general.

Cálculo de variables y razones financieras

library(plm)
## 
## Attaching package: 'plm'
## The following objects are masked from 'package:dplyr':
## 
##     between, lag, lead
#Definimos el dataset uspanel como un data frame:

uspanel <- pdata.frame(uspanel, index= c("firm","q"))
class(uspanel)
## [1] "pdata.frame" "data.frame"

Con el conjunto de datos del panel original (uspanel), creo columnas para las siguientes variables que serán utilizadas para generar las razones financieras más adelante.

Gross profit (grossprofit): Revenue - Cost of good Sold

uspanel$grossprofit = (uspanel$revenue - uspanel$cogs)

Earnings before interest and taxes (ebit): Gross profit - Sales & general administrative expenses

uspanel$ebit = uspanel$grossprofit - uspanel$sgae

Net Income (netincome): ebit - financial expenses - income taxes

uspanel$netincome = uspanel$ebit - uspanel$finexp - uspanel$incometax

Annual market return: ln(adjprice t1) - Ln(adjprice t-4)

uspanel$annual_R = uspanel$adjprice / plm::lag(uspanel$adjprice,k=4) - 1

Market capitalization: (marketcap): original stock price * shareoutstanding.

uspanel$marketcap = uspanel$originalprice * uspanel$sharesoutstanding

Con el mismo conjunto de datos de panel (uspanel), ahora sí, creo columnas para el cálculo de las siguientes razones financieras:

library(psych)
library(statar)

#Cálculo de las razones financieras:

uspanel$roabit = ifelse(plm::lag(uspanel$totalassets,k=4)==0,NA, uspanel$ebit / plm::lag(uspanel$totalassets,k=4))

uspanel$roa = ifelse(plm::lag(uspanel$totalassets,k=4)==0,NA,uspanel$netincome / plm::lag(uspanel$totalassets,k=4))

uspanel$oeps = ifelse(uspanel$sharesoutstanding==0,NA,uspanel$ebit / uspanel$sharesoutstanding)

uspanel$oepsp = ifelse(uspanel$originalprice==0,NA,uspanel$oeps / uspanel$originalprice)

uspanel$bookvalue = uspanel$totalassets - uspanel$totalliabilities

uspanel$bmr = ifelse(uspanel$marketcap==0,NA,uspanel$bookvalue / uspanel$marketcap)

uspanel$annual_r = diff(log(uspanel$adjprice),lag=4)

Estadística descriptiva relevante

nrow(uspanel_2020)
## [1] 1343

En 2020 hubo 1343 empresas públicas activas pertenecientes a la industria de manufactura en el mercado estadounidense.

# Integramos al uspanel_2020 la variable de capitalización de mercado

uspanel_2020$marketcap = uspanel_2020$originalprice * uspanel_2020$sharesoutstanding

# Creamos la variable Top 10 Total Assets para determinar cuáles fueron las empresas más grandes en términos de activos totales de la industria de manufactura en los Estados Unidos. Hacemos lo mismo más adelante para las variables ingreso total y capitalización de mercado.

top_10_tassets <- uspanel_2020 %>%
     arrange(desc(totalassets)) %>%
     top_n(10) %>%
     select(firm, CompanyName, naics1, totalassets, revenue)
## Selecting by marketcap
top_10_tassets

Las empresas más grandes de la industria manufacturera en términos de activos totales fueron Apple, J&J, Pfizer, Intel, P&G, Pepsico, Merck & Co, Coca Cola, Abbott Laboratories y Tesla.

top_10_mcap <- uspanel_2020 %>%
     arrange(desc(marketcap)) %>%
     top_n(10) %>%
     select(firm, CompanyName, naics1, marketcap, revenue)
## Selecting by marketcap
top_10_mcap

Las empresas más grandes de la industria manufacturera en términos de capitalización de mercado fueron Apple, Tesla, J&J, P&G, Coca Cola, Merck & Co, Pepsico, Pfizer, Intel y Abbott Laboratories.

top_10_rev <- uspanel_2020 %>%
     arrange(desc(revenue)) %>%
     top_n(10) %>%
     select(firm, CompanyName, naics1, marketcap, revenue)
## Selecting by marketcap
top_10_rev

Las empresas más grandes de la industria manufacturera en términos de ingresos fueron Apple, Exxon Mobil, GM, Chevron, J&J, GE, Intel, Intl, P&G y Pepsico.

#Analizamos con Q1, Q2 (mediana) y Q3 a la industria en general.

by_manufacturing_summary <- uspanel_2020 %>%
  group_by(naics1) %>%
  summarize(firms = n(),
            median_total_assets = median(totalassets, na.rm = TRUE),
            Q1_total_assets = quantile(totalassets, probs=c(0.25),na.rm = TRUE),
            Q3_total_assets = quantile(totalassets,probs=c(0.75),na.rm=TRUE),
            median_marketcap = median(marketcap, na.rm = TRUE),
            Q1_marketcap = quantile(marketcap, probs=c(0.25),na.rm=TRUE),
            Q3_marketcap = quantile(marketcap,probs=c(0.75),na.rm=TRUE)
            )
by_manufacturing_summary

Hallazgos significativos de la industria manufacturera al año 2020 en el mercado estadounidense:

· El tamaño de empresa típico en la industria manufacturera de las empresas públicas estadounidenses en 2020 en términos de activos totales es de $ 584,954 millones. · El tamaño de empresa típico en la industria manufacturera de las empresas públicas estadounidenses en 2020 en términos de capitalización de mercado es de $ 1,075,660 millones.

· En 2020, aproximadamente el 25% de las empresas públicas pertenecientes a la industria manufacturera de EE. UU. tenían activos totales iguales o inferiores a $ 140,030 millones. · En 2020, aproximadamente el 25% de las empresas públicas pertenecientes a la industria manufacturera de EE. UU. tenían una capitalización de mercado igual o inferior a $ 214,172 millones.

· En 2020, alrededor del 75% de las empresas públicas pertenecientes a la industria manufacturera de EE. UU. tenían activos totales iguales o inferiores a $ 2,804,174 millones. · En 2020, alrededor del 75% de las empresas públicas pertenecientes a la industria manufacturera de EE. UU. tenían una capitalización de mercado igual o inferior a $ 4,378,270 millones.

Visualización de los datos de la industria

#Creamos un histograma para ver más claros los datos y la distribución del market cap en la industria.

firm_selection1 <- uspanel_2020 %>%
            filter(marketcap <=2500000000) %>%
            select(firm, marketcap)
hist(firm_selection1$marketcap, breaks=100)

#Dado que en la industria tenemos a la gigante APPLE, no nos permite visualizar claramente la distribución de las empresas en términos de capitalización de mercado. Por ello, volvemos a elaborar otro histograma recortando a 50000000 el market cap.

firm_selection2 <- uspanel_2020 %>%
            filter(marketcap <=50000000) %>%
            select(firm, marketcap)
hist(firm_selection2$marketcap, breaks=100)

#También podemos visualizar la razón book to market winzorizada de la industria, así como la utilidad operativa por acción deflactada por precio.

uspanel$bmrw <- winsorize(uspanel$bmr,probs=c(0.02,0.98))
## 1.35 % observations replaced at the bottom
## 1.35 % observations replaced at the top
hist(uspanel$bmrw)

uspanel$oepspw <- winsorize(uspanel$oepsp,probs=c(0.02,0.98))
## 1.34 % observations replaced at the bottom
## 1.34 % observations replaced at the top
hist(uspanel$oepspw)

Como podemos ver en este segundo histograma, la capitalización de mercado de las empresas pertenecientes a manufactura en el mercado de los Estados Unidos, mantiene una fuerte y clara oblicuidad a la izquierda. Es decir, que la gran mayoría de las empresas mantiene un capitalización de mercado inferior a 10 mil millones de dólares, billones para los estadounidenses.

Por otra parte, como podemos ver en el gráfico de la utilidad operativa por acción deflactada por precio (oepsp), se mantienen los valores entre un rango de -0.2 a 0.2. A diferencia del rango de las demás industrias (-0.5 a 0.4), el rango de la industria de manufactura es muy bueno, cuando la oepsp es negativo estamos hablando de empresas que arrojan un EBIT negativo. Entonces, en el caso de la industria de manufactura, las empresas que reportaton un -0.2 nos indica que por cada $1.0 invertidos en las acciones, cada accionista podría aportar .20 centavos.

Por último, la razón Book to Market general de la industria de manufactura se mantiene entre un rango de 0 a 0.5; es un excelente resultado para las empresas, ya que ellos justamente buscan estar mejor valuados por el mercado que por sus libros; para así seguir obteniendo financiamiento y expansión.

library(statar)
library(ggplot2)
## 
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
## 
##     %+%, alpha
library(caret)
## Loading required package: lattice
library(Rcpp)

#Winzorizamos la variable marketcap, total assets y revenue dado los extraordinarios valores como APPLE.
#La explicación sobre winzorizar la veremos más adelante.

uspanel_2020$totalassets_w <- winsorize(uspanel_2020$totalassets, probs = c(0,0.95))
## 0.00 % observations replaced at the bottom
## 4.99 % observations replaced at the top
uspanel_2020$revenue_w <- winsorize(uspanel_2020$revenue, probs = c(0,0.95))
## 0.00 % observations replaced at the bottom
## 4.99 % observations replaced at the top
#Creamos una caja de bigotes con la nueva variable winzorizada.

ggplot(uspanel_2020, aes(x=naics1, y=totalassets_w)) +
   geom_boxplot() +
   theme(axis.text.x = element_text(angle = 360, hjust = 1))

ggplot(uspanel_2020, aes(x=naics1, y=revenue_w)) +
   geom_boxplot() +
   theme(axis.text.x = element_text(angle = 360, hjust = 1))

#Podemos volver calcular los histogramas con las variables winzorizadas

hist(uspanel_2020$totalassets_w, breaks=100)

hist(uspanel_2020$revenue_w, breaks=100)

Tanto en las cajas de bigotes como en los histogramas podemos ver una proyección similar, ya que habitualmente la empresas mantienen estructuras verticales muy sólidas y esto propicia a que tanto en capitalización de mercado, activos totales y total de ingresos se mantengan similares y en proporción. Anteriormente vimos el top 10 de empresas en cada una de las variables, sin duda alguna son estas las que nos arroja el programa como los valores máximos observados. El análisis estadístico de estos gráficos lo podemos obtener del apartado de arriba, ya que son los gráficos de las mismas variables.

Modelación con Machine Learning

Machine Learning y tipos de análsis que se pueden hacer con Machine Learning

De acuerdo con BBVA, en su artículo “Machine learning: ¿qué es y cómo funciona?”, el ‘machine learning’ –aprendizaje automático– es una rama de la inteligencia artificial que permite que las máquinas aprendan sin ser expresamente programadas para ello. Una habilidad indispensable para hacer sistemas capaces de identificar patrones entre los datos para hacer predicciones. Esta tecnología está presente en un sinfín de aplicaciones como las recomendaciones de Netflix o Spotify, las respuestas inteligentes de Gmail o el habla de Siri y Alexa. (BBVA, 2021)

En el sentido más básico, Machine Learning (ML) es una forma de implementar inteligencia artificial. Al igual que la IA, el aprendizaje automático es una rama de la informática en la que se diseña o estudia el diseño de algoritmos que pueden aprender. De acuerdo con el sitio oficial de DataCamp, existen varios algoritmos de aprendizaje automático como:

· Decision trees, · Naive Bayes, · Random forest · Support vector machine · K-nearest neighbor, · K-means clustering, · Gaussian mixture model, · Hidden Markov model etc.

Diferencias entre el Machine Learning supervisado y no supervisado

De acuerdo con BiSmart, la modalidad supervisada se podría entender como algoritmos que “aprenden” de los datos introducidos por una persona. En este caso, se necesita la intervención humana para etiquetar, clasificar e introducir los datos en el algoritmo y el algoritmo genera datos de salida esperados, ya que en la entrada han sido etiquetados y clasificados por alguien.

Por su parte, en la modalidad no supervisada no existe la intervención humana en este tipo de Machine Learning. Aquí los algoritmos aprenden de datos con elementos no etiquetados buscando patrones o relaciones entre ellos. En este caso, se introducen datos de entrada sin etiquetar. (BiSmart,2021)

Tipo de Machine Learning utilizado en el análisis y justificación

Para la resolución de este trabajo se han empleado modelos probabilísticos que están ampliamente basados en aplicaciones estadísticas para la evaluación de eventos incontrolables (o factores), así como también la evaluación del riesgo de sus decisiones. Puntualmente empleamos modelos de tipo logit. Los modelos logit, de acuerdo con el Instituto Mora, son un caso particular de los modelos lineales generalizados en los que la distribución es binomial y la función de enlace el logaritmo de las razones de probabilidad, es decir, la función logit que les da nombre. (Instituto More, S.F)

La función glm() de R nos permite ajustar modelos lineales de muchos tipos, incluyendo los que ajustamos con lm(), modelos de Poisson y los logit en los que nos enfocaremos. La sintaxis básica para obtener un modelo lineal es glm(dependiente ~ independiente1 + independiente2, family = binomial(), data=datos").

Los modelos lineales nos permiten predecir el valor de una variable a partir del valor de otra u otras. En su formulación original (Y = mx + B) se desarrollaron para variables dependientes continuas, es decir, para estimar los valores que toma y para cada valor de x a partir de una ordenada al origen y una o más pendientes. No obstante, puede ser de interés tomar como variables dependiente a una categórica y eso es lo que hacen los modelos logit, en su defecto, lo que buscamos en nuestro modelo para este trabajo. Lo que buscaremos predecir son las probabilidades de 0 o 1. Es decir, dado cierto valor de x ¿cuál es la probabilidad de y0 y de y1?

Variables explicativas a utilizar en el análisis

Variable 1: Earnings Per Share (EPS), se calculan las ganancias de una empresa divididas por las acciones en circulación de sus acciones ordinarias. El número resultante sirve como indicador de la rentabilidad de una empresa.

Variabel 2: Book to Market Ratio (BMR), según investopedia, la relación libro-mercado ayuda a los inversores a encontrar el valor de una empresa comparando el valor contable de la empresa con su valor de mercado. Una relación libro-mercado alta podría significar que el mercado está valorando el capital de la empresa a bajo precio en comparación con su valor contable. Muchos inversores están familiarizados con la relación precio / valor contable, que es simplemente la inversa de la fórmula de la relación valor contable / valor mercado. (Investopedia, 2021)

Si el valor de mercado de una empresa cotiza por encima de su valor en libros por acción, se considera que está sobrevalorado. Si el valor en libros es superior al valor de mercado, los analistas consideran que la empresa está infravalorada. La relación libro-mercado se utiliza para comparar el valor liquidativo o valor contable de una empresa con su valor actual o de mercado. Entonces, el nivel de bmr que prefieren todos los ejecutivos es un valor mucho menor que 1.

Si bmr = 0.5 significa que el valor de mercado firme es el doble de su valor en libros. En este caso, la principal diferencia entre el valor de mercado y el valor en libros es la riqueza que la firma ha creado gracias a sus valiosos activos intangibles como el prestigio, la alta calidad y la innovación. De acuerdo con el profesor Alberto Dorantes, la investigación financiera ha descubierto que bmr influye en las ganancias por acción, lo que a su vez influye en los rendimientos actuales y futuros de las acciones. A veces, las empresas con un bmr bajo no experimentan una rentabilidad de las acciones muy alta debido a las altas expectativas de los inversores. Es un hallazgo valioso, pero no universal.

Winsorización de variables realizados

La winsorización es el proceso por el que se sustituye un porcentaje dado de valores extremos en una distribución por lo valores inmediatamente contiguos, con el objetivo de limitar la influencia de estos valores extremos en los estadísticos que resulten de la muestra. En ultima instancia, la winsorización es la transformación de estadísticos limitando los valores extremos en los datos estadísticos para reducir el efecto de valores atípicos, posiblemente falsos.

#Aunque las winsorizamos desde un principio para su visualización, volvemos a traerlas por fines prácticos y de orden

#Probs indicamos la oblicuidad y el percentil al que buscamos acotar nuestros datos.

uspanel$bmrw <- winsorize(uspanel$bmr,probs=c(0.02,0.98))
## 1.35 % observations replaced at the bottom
## 1.35 % observations replaced at the top
uspanel$oepspw <- winsorize(uspanel$oepsp,probs=c(0.02,0.98))
## 1.34 % observations replaced at the bottom
## 1.34 % observations replaced at the top

Manejo de datos para la incorporación del retorno del mercado:

library(quantmod)
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
getSymbols("^GSPC", from="2010-01-01", to= "2020-12-31",
            periodicity="monthly", src="yahoo")
## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
## 
## This message is shown once per session and may be disabled by setting 
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.
## [1] "^GSPC"
QSP500 <- to.quarterly(GSPC,indexAt='startof')
head(QSP500)
##            GSPC.Open GSPC.High GSPC.Low GSPC.Close  GSPC.Volume GSPC.Adjusted
## 2010-01-01   1116.56   1180.69  1044.50    1169.43 279192470000       1169.43
## 2010-04-01   1171.23   1219.80  1028.33    1030.71 354511440000       1030.71
## 2010-07-01   1031.10   1157.16  1010.91    1141.20 260105810000       1141.20
## 2010-10-01   1143.49   1262.60  1131.87    1257.64 257671870000       1257.64
## 2011-01-01   1257.62   1344.07  1249.05    1325.83 240896240000       1325.83
## 2011-04-01   1329.48   1370.58  1258.07    1320.64 245196520000       1320.64
QSP500 = QSP500$GSPC.Adjusted
colnames(QSP500) <- "SP500index"

Ahora podemos calcular el retorno compuesto del mercado

QSP500$SPreturns <- diff(log(QSP500$SP500index),lag=4) 
QSP500.df<-data.frame(q=index(QSP500),coredata(QSP500))

uspanel$q<-as.Date(uspanel$q)

uspanel<-merge(uspanel,QSP500.df,by="q")

Explicar qué variable se utilizó como dependiente y por qué.

La variable dependiente utilizada para el trabajo es el rendimiento anual de las acciones, primero evaluamos el rendimiento contra razones financieras de la empresa para determinar cuáles podrían ser explicativas, y despues nuestra variable dependiente es el rendimiento anual de las acciones pero ahora futuro. Lo que buscamos es determinar con machine learning qué varibles podrían explicar el rendimiento futuro del mercado y así superar el rendimiento promedio histórico del mercado.

Explicar los procesos de: a) entrenamiento del modelo, b) selección de variables explicativas, c) pruebas del modelo entrenado utilizando una muestra de prueba

a) Selección de variables explicativas

Primero que nada, ejecutamos un modelo de regresión múltiple para examinar si la utilidad operativa por acción deflactada por precio (oepsp) y la rentabilidad operativa sobre los activos (roabit) están relacionados con los rendimientos anuales contemporáneos. No comprobamos todavía si oepsp y roabit están relacionados con los rendimientos un trimestre después; solo comprobamos si están relacionados con los rendimientos en el mismo período para todos los trimestres de empresas.

uspanelactive <- uspanel[uspanel$status == "active",]

#lm = linear model

LM1 <- lm(annual_r ~ oepsp + roabit, data = uspanelactive, na.action = na.omit)
summary(LM1)
## 
## Call:
## lm(formula = annual_r ~ oepsp + roabit, data = uspanelactive, 
##     na.action = na.omit)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.3513 -0.2358  0.0366  0.2668 10.8563 
## 
## Coefficients:
##               Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)  0.0050164  0.0027186   1.845              0.0650 .  
## oepsp        0.6751529  0.0079798  84.608 <0.0000000000000002 ***
## roabit      -0.0005494  0.0002213  -2.482              0.0131 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5301 on 38563 degrees of freedom
##   (25982 observations deleted due to missingness)
## Multiple R-squared:  0.1566, Adjusted R-squared:  0.1565 
## F-statistic:  3580 on 2 and 38563 DF,  p-value: < 0.00000000000000022

Interpretación:

El siguiente modelo de regresión lineal contempla como variable dependiente el retorno anual de las acciones (trimestre a trimestre) y como variables independientes la utilidad operativa por acción deflactada por precio (oepsp) y la rentabilidad operativa sobre los activos (roabit).

El coeficiente beta de oepsp es de 0.6751529. Este coeficiente es positivo y estadísticamente significativo ya que su valor p (0.0000000000000002) es mucho menor que 0.05. Esto nos indica que al menos el 95% de las veces, en un futuro próximo este coeficiente será positivo. En otras palabras, existe una relación significativa y positiva entre las ganancias operativas por acción deflactadas por el precio y los rendimientos anuales de las acciones después de considerar el efecto del rendimiento operativo sobre los activos sobre los rendimientos anuales de las acciones. Es decir, después de considerar el efecto de roabit, por cada cambio de +1 en oepsp, se espera que el cambio en el rendimiento anual sea de aproximadamente 0.6751529.

Por otra parte, el coeficiente de roa operacional (roabit) es de -0.0005494; este coeficiente es negativo y estadísticamente significativo ya que su valor p (0.0131) es menor que 0.05. Por lo tanto, después de considerar el efecto de oeps en el rendimiento de las acciones, podemos determinar que sí hay un efecto relevante o relación entre el roabit y la rentabilidad anual de las acciones. Cabe resaltar que a diferencia de otras industrias, en la manufactura se logra apreciar esta relación, en algunas otras no hay soporte estadístico para concluir lo mismo.

Conluimos que en efecto, si hay una relación significativa entre la utilidad operativa por acción deflactada por precio (oepsp) y la rentabilidad operativa sobre los activos (roabit) con relación a los rendimientos anuales. Ahora bien, hacemos uso de la variable dependiente de esta regresión múltiple para que sea el rendimiento anual de las acciones un trimestre futuro.

La función plm nos permite usar los valores futuros de una variable en los datos del panel:

PLM1 <- plm(plm::lag(annual_r,-1) ~ oepsp + roabit, data=uspanelactive, model="pooling", na.action = na.omit)
summary(PLM1)
## Pooling Model
## 
## Call:
## plm(formula = plm::lag(annual_r, -1) ~ oepsp + roabit, data = uspanelactive, 
##     na.action = na.omit, model = "pooling")
## 
## Unbalanced Panel: n = 40, T = 61-1194, N = 29238
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -4.456315 -0.221499  0.074515  0.303938  3.633940 
## 
## Coefficients:
##                  Estimate    Std. Error t-value              Pr(>|t|)    
## (Intercept) -0.0323972521  0.0034150712 -9.4866 < 0.00000000000000022 ***
## oepsp        0.0416315419  0.0095916346  4.3404            0.00001427 ***
## roabit      -0.0000021233  0.0000030303 -0.7007                0.4835    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    9822.3
## Residual Sum of Squares: 9815.9
## R-Squared:      0.00066072
## Adj. R-Squared: 0.00059235
## F-statistic: 9.66446 on 2 and 29235 DF, p-value: 0.000063704

Interpretación:

Para este segundo modelo, el coeficiente beta de oepsp es de 0.0416315419. Este coeficiente es positivo y estadísticamente significativo ya que su valor p (0.0000000000000002) es mucho menor que 0.05. Esto nos indica que al menos el 95% de las veces, en un futuro próximo este coeficiente será positivo. En otras palabras, existe una relación significativa y positiva entre las ganancias operativas por acción deflactadas por el precio y los rendimientos anuales de las acciones después de considerar el efecto del rendimiento operativo sobre los activos sobre los rendimientos anuales de las acciones. Es decir, después de considerar el efecto de roabit, por cada cambio de +1 en oepsp, se espera que el cambio en el rendimiento anual sea de aproximadamente 0.0416315419.

Por otra parte, el coeficiente de roa operacional (roabit) es de -0.000002123; este coeficiente es negativo y estadísticamente no significativo ya que su valor p (0.4835) es mayor que 0.05. Por lo tanto, después de considerar el efecto de oeps en el rendimiento de las acciones, podemos determinar que sí hay un efecto relevante o relación entre el roabit y la rentabilidad anual de las acciones.

Conluimos que en efecto, si hay una relación significativa entre la utilidad operativa por acción deflactada por precio (oepsp) y la proyección de los rendimiento de las acciones, sin embargo, la rentabilidad operativa sobre los activos (roabit) no proyecta una relación significativa con los rendimientos anuales.

# plm = panel linear model
# La función plm nos permite usar los valores rezagados o futuros de una variable en los datos del panel.
# Variables winsorizadas.

uspanelactive$size = log(uspanelactive$marketcap)

PLM2 <- plm(plm::lag(annual_r,-1) ~ bmrw + oepspw + size, data=uspanelactive, model="pooling", na.action = na.omit)
summary(PLM2)
## Pooling Model
## 
## Call:
## plm(formula = plm::lag(annual_r, -1) ~ bmrw + oepspw + size, 
##     data = uspanelactive, na.action = na.omit, model = "pooling")
## 
## Unbalanced Panel: n = 40, T = 61-1232, N = 30971
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -5.050723 -0.221188  0.073505  0.301620  3.616739 
## 
## Coefficients:
##                Estimate  Std. Error t-value      Pr(>|t|)    
## (Intercept) -0.00309411  0.02626108 -0.1178        0.9062    
## bmrw        -0.05065529  0.00881081 -5.7492 0.00000000905 ***
## oepspw       0.12924974  0.02449852  5.2758 0.00000013305 ***
## size        -0.00036918  0.00178967 -0.2063        0.8366    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    10362
## Residual Sum of Squares: 10335
## R-Squared:      0.0025933
## Adj. R-Squared: 0.0024966
## F-statistic: 26.8382 on 3 and 30967 DF, p-value: < 0.000000000000000222

Tras haber determinado que el oepsp si mantiene un relación significativa con los rendimientos futuros de la acciones y que existe soporte estadístico que avale esta variable; la incluímos a nuestro tercer modelo, y en esta tercer ocasión agregamos la razón financiero BMR. Cabe destacar, que por fines de gestión de data se decidió winsorizar ambas variables y de este modo validar estadísticamente su uso.

En este tercer modelo, obtenemos que tanto el BMR como el OEPSP mantienen una relación significativa y soportada estadísticamente dado su valor p. Sin embargo, el crecimiento de la empresa no resulta ser significativo.

b) Entrenamiento del modelo

Ya que tenemos una idea más clara sobre qué variables podrían explicar nuestros posibles rendimientos futuros, es momento de desarrollar nuestro modelo final para vencer al mercado.

#Variable rendimiento superior al del S&P 500
uspanelactive$r_AboveMarket = ifelse(uspanelactive$annual_r>uspanelactive$SPreturns, 1, 0)

uspanelactive$F4_rAboveMarket = plm::lag(uspanelactive$r_AboveMarket, -4)

table(uspanelactive$r_AboveMarket)
## 
##     0     1 
## 22875 17685

Número de periodos en el que las empresas de manufactura superaron el rendimiento del mercado (1).

#Modelo final, variables winsorizadas, proyección a vencer el mercado.
# Size = log(MarketCap)

uspanelactive$size = log(uspanelactive$marketcap)

GLM1 <- glm(F4_rAboveMarket ~ bmrw + oepspw + size, data=uspanelactive, family = "binomial", na.action = na.omit)
summary(GLM1)
## 
## Call:
## glm(formula = F4_rAboveMarket ~ bmrw + oepspw + size, family = "binomial", 
##     data = uspanelactive, na.action = na.omit)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9287  -1.0743  -0.4518   1.0870   3.3489  
## 
## Coefficients:
##              Estimate Std. Error z value             Pr(>|z|)    
## (Intercept) -0.266935   0.089746  -2.974              0.00294 ** 
## bmrw        -1.683591   0.037769 -44.575 < 0.0000000000000002 ***
## oepspw       4.174284   0.126051  33.116 < 0.0000000000000002 ***
## size         0.055810   0.006009   9.288 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 53938  on 39344  degrees of freedom
## Residual deviance: 48251  on 39341  degrees of freedom
##   (25203 observations deleted due to missingness)
## AIC: 48259
## 
## Number of Fisher Scoring iterations: 5

Todas las variables independientes mantienen una relación significativa y soportada estadísitcamente a los futuros rendimientos del mercado dado su valor p inferior a 0.05.

#Creación nuevo panel para TRAINING Y TESTING.

library(caret)
library(Rcpp)

uspanel_proy <- uspanelactive %>%
  select(firm, q, year, fiscalmonth, annual_r, F4_rAboveMarket, bmrw, oepspw, size) %>%
  na.omit() %>%
  as.data.frame()

c) Pruebas del modelo entrenado utilizando una muestra de prueba

set.seed(1221)

#TRAINING

trainedata = createDataPartition(uspanel_proy$F4_rAboveMarket, p=.80, list=FALSE)

training = uspanel_proy[trainedata,]
testing = uspanel_proy[-trainedata,]

#Training del modelo final con variables explicativas definitivas.

GLM_ml <- glm(F4_rAboveMarket ~ bmrw + oepspw + size, data=training, family = "binomial", na.action = na.omit)
summary(GLM_ml)
## 
## Call:
## glm(formula = F4_rAboveMarket ~ bmrw + oepspw + size, family = "binomial", 
##     data = training, na.action = na.omit)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9086  -1.0763  -0.4659   1.0925   3.3293  
## 
## Coefficients:
##              Estimate Std. Error z value             Pr(>|z|)    
## (Intercept) -0.219428   0.100466  -2.184                0.029 *  
## bmrw        -1.662595   0.042277 -39.326 < 0.0000000000000002 ***
## oepspw       4.155218   0.140796  29.512 < 0.0000000000000002 ***
## size         0.051103   0.006715   7.610   0.0000000000000275 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 43137  on 31475  degrees of freedom
## Residual deviance: 38739  on 31472  degrees of freedom
## AIC: 38747
## 
## Number of Fisher Scoring iterations: 5

Al haber obtenido estos resultados, podemos observar que las tres variables empleadas (BMR, OEPSP y SIZE) mantienen una relación significativa a la proyección de los retornos de las acciones. Por lo que comenzaremos con la fase TESTING.

#Evaluamos el nivel de precisión del modelo con los casos que no se incluyeron en la calibración.
#Evaluamos qué tan bien predice nuestro modelo final.

testing <- testing %>%
  mutate(pred=predict.glm(GLM_ml,newdata=testing,type=c("response")))

#Probabilidad de vencer el retorno del mercado con el modelo final.

hist(testing$pred)

La probabilidad de vencer al mercado con este modelo es de 45% - 65%, excelente porcentaje.

#TESTING

#Indicamos probabilidad superior a 60%

testing$F4_rAboveMarket_prediction =ifelse(testing$pred > 0.5,1,0)
table(testing$F4_rAboveMarket_prediction)
## 
##    0    1 
## 4467 3402
#Generamos variables factor para la confusionMatrix

testing$F4_rAboveMarket_pred_f_ <- factor(testing$F4_rAboveMarket)
testing$F4_rAboveMarket_pred_f_pred <- factor(testing$F4_rAboveMarket_prediction)

#Predicción 1 y 0 vs. si realmente le ganó.

confusionMatrix(testing$F4_rAboveMarket_pred_f_pred, testing$F4_rAboveMarket_pred_f_, positive="1")
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 3104 1363
##          1 1292 2110
##                                              
##                Accuracy : 0.6626             
##                  95% CI : (0.652, 0.673)     
##     No Information Rate : 0.5586             
##     P-Value [Acc > NIR] : <0.0000000000000002
##                                              
##                   Kappa : 0.3143             
##                                              
##  Mcnemar's Test P-Value : 0.1743             
##                                              
##             Sensitivity : 0.6075             
##             Specificity : 0.7061             
##          Pos Pred Value : 0.6202             
##          Neg Pred Value : 0.6949             
##              Prevalence : 0.4414             
##          Detection Rate : 0.2681             
##    Detection Prevalence : 0.4323             
##       Balanced Accuracy : 0.6568             
##                                              
##        'Positive' Class : 1                  
##

El nivel de precisión de nuestro modelo es de 65.32%, la predicción del modelo y la realidad coinciden en 2103 veces. Es decir, coincide entre la predicción del modelo y el resultado real. Para la complejidad que conlleva crear un modelo que logre predecir el futuro rendimiento de algun activo, es muy buen resultado.

#Modelo final elaborado con las variables explicativas BMRW OEPSPW y SIZE.

ModeloFinal <- glm(F4_rAboveMarket ~ bmrw + oepspw + size, data=uspanel_proy, family = "binomial", na.action = na.omit)
summary(ModeloFinal)
## 
## Call:
## glm(formula = F4_rAboveMarket ~ bmrw + oepspw + size, family = "binomial", 
##     data = uspanel_proy, na.action = na.omit)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9287  -1.0743  -0.4518   1.0870   3.3489  
## 
## Coefficients:
##              Estimate Std. Error z value             Pr(>|z|)    
## (Intercept) -0.266935   0.089746  -2.974              0.00294 ** 
## bmrw        -1.683591   0.037769 -44.575 < 0.0000000000000002 ***
## oepspw       4.174284   0.126051  33.116 < 0.0000000000000002 ***
## size         0.055810   0.006009   9.288 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 53938  on 39344  degrees of freedom
## Residual deviance: 48251  on 39341  degrees of freedom
## AIC: 48259
## 
## Number of Fisher Scoring iterations: 5

Una vez más observamos como todas las variables son soportadas estadísticamente hablando y significativas contra el retorno futuro del mercado.

library(dplyr)

#Simplificamos nuestro panel para 

uspanel_proy <- pdata.frame(uspanel_proy, index = c("firm","q"))
lastq <- uspanel_proy %>%
      filter(q=="2020-10-01") %>%
      select(firm, annual_r, F4_rAboveMarket, bmrw, oepspw, size)
#Probabilidad de vencer al mercado con nuestra modelación

#Creamos el modelo de predicción con GLM

lastq <- lastq %>%
  mutate(pred=predict.glm(ModeloFinal, newdata=lastq, type=c("response")))

hist(lastq$pred)

Probabilidad del modelo cerca del 60% para vencer al mercado.

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ tibble  3.1.0     ✓ purrr   0.3.4
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x ggplot2::%+%()   masks psych::%+%()
## x ggplot2::alpha() masks psych::alpha()
## x plm::between()   masks dplyr::between()
## x dplyr::filter()  masks stats::filter()
## x xts::first()     masks dplyr::first()
## x plm::lag()       masks dplyr::lag(), stats::lag()
## x xts::last()      masks dplyr::last()
## x plm::lead()      masks dplyr::lead()
## x purrr::lift()    masks caret::lift()
#Seleccionamos el top 10 de nuestras predicciones que vencen al mercado para el futuro con base a la predicción del 4Q y en orden descendente.

Top10_Man_Firms = lastq%>%
  arrange(desc(lastq$pred)) %>%
  top_n(10)
## Selecting by pred
Top10_Man_Firms

Estas son las firmas que nuestro modelo predice que venceran al mercado, es decir, nuestro portafolio ganador en la industria de manufactura. Ahora elaboraremos el portafolio y evaluaremos la proyección del rendimiento.

#Vectorizamos nuestro top 10
Portafolio = as.vector(Top10_Man_Firms$firm)
Portafolio
##  [1] "VGR" "EAF" "RRD" "HPQ" "TUP" "PM"  "MO"  "SUP" "PBI" "BW"
#Traemos sus valores de Yahoo Finance para determinar su retorno
getSymbols(Symbols=Portafolio, from="2021-01-01", to="2021-10-15")
## pausing 1 second between requests for more than 5 symbols
## pausing 1 second between requests for more than 5 symbols
## pausing 1 second between requests for more than 5 symbols
## pausing 1 second between requests for more than 5 symbols
## pausing 1 second between requests for more than 5 symbols
## pausing 1 second between requests for more than 5 symbols
##  [1] "VGR" "EAF" "RRD" "HPQ" "TUP" "PM"  "MO"  "SUP" "PBI" "BW"
Portafolio_m = merge(VGR, EAF, RRD, HPQ, TUP, PM, MO, SUP, PBI, BW)

#Extraemos su precio ajustado
Portafolio_m = Ad(Portafolio_m)

#Calculamos el retorno de las empresas del Top10
r_calculado = (Portafolio_m/stats::lag(Portafolio_m))-1
r_calculado = na.omit(r_calculado)


rE_Portafolio = apply(r_calculado, MARGIN=1, FUN = mean)

#utilizamos este paquete que cuenta con la función de como se vería $1 con cierto rendimiento.
library(PerformanceAnalytics)

charts.PerformanceSummary(rE_Portafolio,
                          main = "Rendimiento del portafolio",
                          wealth.index = TRUE)

# Análisis de resultado, recomendación final y conclusión

Tras haber elaborado el modelo, obtenemos nuestro resultados finales. En esta ocasión el portafolio ganador como vimos incluye las acciones de VGR, EAF, RRD, HPQ, TUP, PM, MO, SUP, PBI y BW. Nuestro modelo proyecta un rendimiento aproximado el 45% anual. Cabe destacar que el portafolio no se encuentra optimizado y se asume un peso igual sobre todos los activos.

En mi opinión, es un rendimiento extraordinario, sin embargo, esto no es más que un modelo y se emplearon predicciones, es decir, no es algo definitivo. Como bien sabemos rendimientos históricos no garantizan rendimientos futuros. Por otra parte, me gustaría concluir que yo esperaba la selección de activos como Apple, Tesla, Intel, etc. ya que como vimos en un principio son las empresas más grandes de la industria. Habitualmente son empresas que vienen venciendo al mercado y por alguna razón, quizás un error en el modelo, no decidió elegir estas acciones.

#Referencias:

https://www.investopedia.com/terms/e/eps.asp

https://www.investopedia.com/terms/b/booktomarketratio.asp

https://blog.bismart.com/es/diferencias-machine-learning-supervisado-no-supervisado

https://www.institutomora.edu.mx/testu/sitepages/martinpaladino/modelos_logit_con_r.html