Code for Session 8

Note that the directories used to store data are likely different on your computer, and such references will need to be changed before using any such code.

library(knitr)
library(kableExtra)
html_df <- function(text, cols=NULL, col1=FALSE, full=F) {
  if(!length(cols)) {
    cols=colnames(text)
  }
  if(!col1) {
    kable(text,"html", col.names = cols, align = c("l",rep('c',length(cols)-1))) %>%
      kable_styling(bootstrap_options = c("striped","hover"), full_width=full)
  } else {
    kable(text,"html", col.names = cols, align = c("l",rep('c',length(cols)-1))) %>%
      kable_styling(bootstrap_options = c("striped","hover"), full_width=full) %>%
      column_spec(1,bold=T)
  }
}

library(tidyverse)

readRDS("../../Data/corp_summary.rds")

corp_FOG <- readRDS("../../Data/corp_readability_FOG.rds")
corp_FOG %>%
  head() %>%
  html_df()

document	FOG
0000002178-14-000010.txt	21.03917
0000003499-14-000005.txt	20.36549
0000003570-14-000031.txt	22.24386
0000004187-14-000020.txt	18.75720
0000004457-14-000036.txt	19.22683
0000004904-14-000019.txt	20.51594

summary(corp_FOG$FOG)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  14.33   20.32   21.01   21.05   21.75   35.37 

ggplot(corp_FOG, aes(x=FOG)) + geom_density()

df_SIC <- read.csv('../../Data/Filings2014.csv') %>%
  select(accession, regsic) %>%
  mutate(accession=paste0(accession, ".txt")) %>%
  rename(document=accession) %>%
  mutate(industry = case_when(
    regsic >=0100 & regsic <= 0999 ~ "Agriculture",
    regsic >=1000 & regsic <= 1499 ~ "Mining",
    regsic >=1500 & regsic <= 1799 ~ "Construction",
    regsic >=2000 & regsic <= 3999 ~ "Manufacturing",
    regsic >=4000 & regsic <= 4999 ~ "Utilities",
    regsic >=5000 & regsic <= 5199 ~ "Wholesale Trade",
    regsic >=5200 & regsic <= 5999 ~ "Retail Trade",
    regsic >=6000 & regsic <= 6799 ~ "Finance",
    regsic >=7000 & regsic <= 8999 ~ "Services",
    regsic >=9100 & regsic <= 9999 ~ "Public Admin" )) %>%
  group_by(document) %>%
  slice(1) %>%
  ungroup()
corp_FOG <- corp_FOG %>% left_join(df_SIC)

Joining, by = "document"

corp_FOG %>%
  head() %>%
  html_df()

document	FOG	regsic	industry
0000002178-14-000010.txt	21.03917	5172	Wholesale Trade
0000003499-14-000005.txt	20.36549	6798	Finance
0000003570-14-000031.txt	22.24386	4924	Utilities
0000004187-14-000020.txt	18.75720	4950	Utilities
0000004457-14-000036.txt	19.22683	7510	Services
0000004904-14-000019.txt	20.51594	4911	Utilities

ggplot(corp_FOG[!is.na(corp_FOG$industry),], aes(x=factor(industry), y=FOG)) +
  geom_violin(draw_quantiles = c(0.25, 0.5, 0.75)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplot(corp_FOG[!is.na(corp_FOG$industry),], aes(x=FOG)) +
  geom_density() + facet_wrap(~industry)

library(lattice)
densityplot(~FOG | industry,
            data=corp_FOG,
            plot.points=F,
            main="Fog index distibution by industry (SIC)",
            xlab="Fog index",
            layout=c(3,3))

library(DT)
readRDS('../../Data/corp_kwic.rds') %>%
  mutate(text=paste(pre,keyword,post)) %>%
  left_join(select(df_SIC, document, industry), by = c("docname" = "document")) %>%
  select(docname, industry, text) %>%
  datatable(options = list(pageLength = 5), rownames=F)

gw_count <- readRDS('../../Data/corp_kwic.rds') %>%
  left_join(select(df_SIC, document, industry), by = c("docname" = "document")) %>%
  group_by(industry) %>%
  mutate(docs_gw = n()) %>%
  slice(1) %>%
  ungroup() %>%
  select(industry, docs_gw)

corp_FOG %>% group_by(industry) %>%
  mutate(docs = n()) %>%
  slice(1) %>%
  ungroup() %>%
  left_join(gw_count) %>%
  mutate(docs_gw = ifelse(is.na(docs_gw),0,docs_gw)) %>%
  mutate(`Global warming` = docs_gw / docs,
         `Does not mention` = (docs - docs_gw) / docs) %>%
  gather(mention, percent, `Global warming`, `Does not mention`) %>%
  ggplot() + 
  geom_col(aes(x = industry, y = percent, fill=mention)) + 
  ylab("Percent of annual reports") +
  xlab("Industry (SIC code)") +
  ggtitle("Industries discussing global warming in 2014") + 
  scale_y_continuous(labels = scales::percent) + 
  scale_fill_manual(values=c("#CCCCCCAA", "#88CC88")) +
  theme(axis.text.x = element_text(angle=45, hjust=1))

readRDS('../../Data/corp_dfm_feat.rds')

$`Wholesale Trade`
compani    oper million financi product 
  30371   20340   18085   17552   17300 

$Finance
compani    loan financi  decemb million 
 438185  392164  299978  286791  274376 

$Utilities
   oper million compani financi  includ 
 112038  107322  101971   79010   76604 

$Services
compani    oper million financi  servic 
 222276  145506  138397  131881  120817 

$Manufacturing
compani product million    oper financi 
 434805  368900  275829  240181  231687 

$Mining
compani    oper     gas     oil  decemb 
  97798   92076   74150   65532   60475 

$Construction
compani million    oper financi  decemb 
  15479   14885   12431   10899   10149 

$`Retail Trade`
compani    oper million financi  includ 
  62780   43637   41428   35824   32478 

$Agriculture
compani    oper financi    year product 
   4200    3016    2949    2756    2750 

$<NA>
numeric(0)

readRDS('../../Data/corp_tfidf_feat.rds')

$`Wholesale Trade`
  graybar  grainger       oil   million     bottl 
0.3140485 0.2899255 0.2187512 0.2184815 0.2122642 

$Finance
       ab   mortgag depositor      loan      reit 
 9.863862  7.414096  6.192815  5.109854  5.046502 

$Utilities
     gas      fcc  pipelin   energi aircraft 
2.005220 1.484092 1.227766 1.164767 1.020255 

$Services
    game   client   casino  million  softwar 
2.394468 1.760647 1.635549 1.496073 1.404740 

$Manufacturing
  clinic      fda    trial     drug  patient 
7.057913 5.487707 3.949705 3.935010 3.799611 

$Mining
     gas      oil    drill     well   explor 
6.550322 6.308205 4.935983 2.412994 2.035304 

$Construction
homebuild      home     iveda      layn       alp 
0.5143533 0.3827212 0.3557692 0.2360279 0.2303252 

$`Retail Trade`
   restaur      store merchandis  franchise   franchis 
 2.6829714  1.5131383  1.3382872  0.8695339  0.5706876 

$Agriculture
      yew      uspb  mushroom       prc       nbp 
0.2894138 0.2140732 0.2066838 0.2031097 0.1862960 

$<NA>
numeric(0)

readRDS('../../Data/corp_tfidf_bank.rds')

         ab     mortgag   depositor        loan        reit       trust 
   9.863862    7.414096    6.192815    5.109854    5.046502    4.394811 
    reinsur      truste       estat      tenant    instruct partnership 
   3.809024    3.607591    3.188824    3.100092    2.970419    2.697215 
       real     million        pool        fdic   residenti     bancorp 
   2.506670    2.482285    2.287610    2.238533    2.149133    2.074819 
    obligor        rmbs 
   2.055811    2.055453 

out <- readRDS('../../Data/corp_out_stm.rds')
#out <- readRDS('../../Data/corp_out_lda.rds')
#out <- readRDS('../../Data/corp_out_topicmodels.rds')

out$documents[[1]][,386:390]

      [,1]  [,2]  [,3]  [,4]  [,5]
[1,] 14590 14593 14598 14614 14625
[2,]     1     1    38     3     1

out$vocab[c(out$documents[[1]][,386:390][1,])]

[1] "earlier"  "earliest" "earn"     "earthen"  "eas"     

library(stm)
topics <- readRDS('../../Data/corp_stm_topics.rds')

labelTopics(topics)

Topic 1 Top Words:
     Highest Prob: properti, oper, million, decemb, compani, interest, leas 
     FREX: ffo, efih, efh, tenant, hotel, casino, guc 
     Lift: aliansc, baluma, change-of-ownership, crj700s, directly-reimburs, escena, hhmk 
     Score: reit, hotel, game, ffo, tenant, casino, efih 
Topic 2 Top Words:
     Highest Prob: compani, stock, share, common, financi, director, offic 
     FREX: prc, asher, shaanxi, wfoe, eit, hubei, yew 
     Lift: aagc, abramowitz, accello, akash, alix, alkam, almati 
     Score: prc, compani, penni, stock, share, rmb, director 
Topic 3 Top Words:
     Highest Prob: product, develop, compani, clinic, market, includ, approv 
     FREX: dose, preclin, nda, vaccin, oncolog, anda, fdas 
     Lift: 1064nm, 12-001hr, 25-gaug, 2ml, 3shape, 503b, 600mg 
     Score: clinic, fda, preclin, dose, patent, nda, product 
Topic 4 Top Words:
     Highest Prob: invest, fund, manag, market, asset, trade, interest 
     FREX: uscf, nfa, unl, uga, mlai, bno, dno 
     Lift: a-1t, aion, apx-endex, bessey, bolduc, broyhil, buran 
     Score: uscf, fhlbank, rmbs, uga, invest, mlai, ung 
Topic 5 Top Words:
     Highest Prob: servic, report, file, program, provid, network, requir 
     FREX: echostar, fcc, fccs, telesat, ilec, starz, retransmiss 
     Lift: 1100-n, 2-usb, 2011-c1, 2012-ccre4, 2013-c9, aastra, accreditor 
     Score: entergi, fcc, echostar, wireless, broadcast, video, cabl 
Topic 6 Top Words:
     Highest Prob: loan, bank, compani, financi, decemb, million, interest 
     FREX: nonaccru, oreo, tdrs, bancorp, fdic, charge-off, alll 
     Lift: 100bp, 4-famili, acnb, acquired-impair, amerihom, ameriserv, annb 
     Score: fhlb, loan, bank, mortgag, risk-weight, tdrs, nonaccru 
Topic 7 Top Words:
     Highest Prob: compani, million, oper, financi, revenu, result, includ 
     FREX: vmware, imax, franchise, merchandis, affinion, exhibitor, softwar 
     Lift: 4.75x, 9corpor, accessdm, acvc, adtech, adxpos, aecsoft 
     Score: million, product, restaur, custom, game, video, merchandis 
Topic 8 Top Words:
     Highest Prob: compani, insur, million, loss, financi, invest, rate 
     FREX: policyhold, reinsur, lae, dac, annuiti, ambac, cede 
     Lift: agcpl, ahccc, amcareco, argoglob, asil, connecticut-domicil, feinsod 
     Score: reinsur, policyhold, lae, onebeacon, insur, million, dac 
Topic 9 Top Words:
     Highest Prob: million, compani, oper, financi, cost, product, decemb 
     FREX: wafer, alcoa, pepco, dpl, nstar, usec, kcsm 
     Lift: 1.5mw, 11-hour, 1ynanomet, 3dfx, 3ms, 3pd, 40g 
     Score: million, product, ameren, cleco, manufactur, wafer, postretir 
Topic 10 Top Words:
     Highest Prob: gas, oper, oil, natur, million, cost, decemb 
     FREX: ngl, ngls, oneok, mgp, permian, qep, wes 
     Lift: 12asset, 1businesscommerci, 94-mile, aivh, amargo, amopp, angell 
     Score: gas, drill, oil, ngl, crude, unithold, ngls 

out$meta$industry <- factor(out$meta$industry)

doc_topics = data.frame(document=names(out$documents),
                        industry=out$meta$industry,
                        topic=1,
                        weight=topics$theta[,1])
for (i in 2:10) {
  temp = data.frame(document=names(out$documents),
                    industry=out$meta$industry,
                    topic=i,
                    weight=topics$theta[,i])
  doc_topics = rbind(doc_topics, temp)
}

# Proportional topics (%)
doc_topics <- doc_topics %>%
  group_by(document) %>%
  mutate(topic_prop = weight / sum(weight)) %>%
  ungroup()

# Manually label topics
topic_labels = data.frame(topic = 1:10,
                          topic_name = c('Real Estate', 'Management', 'Product',
                                         'Investment', 'Services', 'Financing',
                                         'Service2', 'Insurance', 'Industrial',
                                         'Utility'))

doc_topics <- doc_topics %>% left_join(topic_labels)

doc_topics %>% filter(document=='0001104659-14-015152.txt')

doc_topics %>%
  filter(document=='0001104659-14-015152.txt' |
         document=='0000019617-14-000289.txt') %>%
  mutate(Company=ifelse(document=='0001104659-14-015152.txt', 'Citi','JPM')) %>%
  ggplot(aes(x=factor(topic_name), y=topic_prop, fill=factor(topic_name))) + 
  geom_col() + facet_wrap(~Company) + 
  theme(axis.text.x=element_blank(),axis.ticks.x = element_blank())

doc_topics %>%
  group_by(industry, topic) %>%
  mutate(topic_prop = mean(topic_prop)) %>%
  slice(1) %>%
  ungroup() %>%
  ggplot(aes(x=factor(topic_name), y=topic_prop, fill=factor(topic_name))) + 
  geom_col() + facet_wrap(~industry) +
  theme(axis.text.x=element_blank(),axis.ticks.x = element_blank())

library(tidyr)
wide_topics <- spread(doc_topics[,c(1,2,5,6)], topic_name, topic_prop)
mat <- wide_topics[,3:12]

mat[,1:6] %>% head() %>% html_df()

Financing	Industrial	Insurance	Investment	Management	Product
0.0105862	0.1578543	0.1088631	0.0004632	0.1161191	0.0002101
0.0467173	0.0059438	0.0235389	0.0005284	0.0801189	0.0001432
0.0069105	0.0351987	0.0003661	0.0201215	0.0023672	0.0000186
0.0870371	0.8271759	0.0003259	0.0003334	0.0206444	0.0000485
0.0036086	0.2680866	0.2677154	0.0008808	0.0026448	0.0000949
0.0000976	0.5299432	0.0001593	0.0007533	0.0009532	0.0000318

set.seed(6845868)
clusters <- kmeans(mat, 9)
clusters$cluster %>% head()

[1] 7 3 2 9 4 7

cbind(as.data.frame(clusters$center), data.frame(kmean=1:9)) %>%
  gather("Topics","weights",-kmean) %>%
  ggplot(aes(x=factor(Topics), y=weights, fill=factor(Topics))) +
  geom_col() + 
  facet_wrap(~kmean) + 
  theme(axis.text.x=element_blank(),axis.ticks.x = element_blank())

library(cluster) # Uses PCA (principle component analysis)
clusplot(mat, clusters$cluster, color=TRUE, shade=TRUE, 
         labels=4)

library(Rtsne)
dups <- which(duplicated(mat))
wide_nodup <- wide_topics[-dups,]
wide_nodup$kmean <- clusters$cluster[-dups]
tsne_data <- readRDS('../../Data/corp_tsne.rds')

wide_nodup <- wide_nodup %>%
  mutate(tsne1 = tsne_data$Y[, 1], tsne2 = tsne_data$Y[, 2])

ggplot(wide_nodup, aes(x = tsne1, y = tsne2, colour = industry)) + 
    geom_point(alpha = 0.3) + theme_bw()

ggplot(wide_nodup, aes(x = tsne1, y = tsne2, colour = factor(kmean))) + 
    geom_point(alpha = 0.3) + theme_bw()

ggplot(wide_nodup, aes(x=kmean)) + geom_bar() + facet_wrap(~factor(industry))

ggplot(wide_nodup, aes(x=tsne1, y=tsne2, color=factor(kmean))) + geom_point() +
  facet_wrap(~factor(industry))

ggplot(wide_nodup, aes(x=tsne1, y=tsne2, color=factor(industry))) + geom_point() +
  facet_wrap(~factor(kmean))

#wide_topics$dist <- sqrt(rowSums(mat - fitted(clusters))^2)
wide_topics$dist <- sqrt(rowSums(abs(mat - fitted(clusters))))
wide_topics[,c(1,2,3,5,13)] %>% arrange(desc(dist)) %>% slice(1:5) %>% html_df()

document	industry	Financing	Insurance	dist
0001171500-14-000007.txt	Finance	0.0003177	0.9972499	1.341244
0001193125-14-077320.txt	Finance	0.0001725	0.9794704	1.337283
0001095073-14-000008.txt	Finance	0.0002339	0.9916079	1.337031
0000356130-14-000052.txt	Finance	0.0002991	0.9845263	1.334780
0000021175-14-000021.txt	Finance	0.0036298	0.9875105	1.333963

wide_topics[,c(2,5,8,9,10,11,13)] %>% filter(industry!="Finance") %>% arrange(desc(dist)) %>% mutate(id=1:n())%>% select(id,everything()) %>% slice(1:2,7,6,8) %>% html_df()

id	industry	Insurance	Product	Real Estate	Service2	Services	dist
1	Services	0.5161854	0.2641739	0.1112599	0.0136804	0.0764400	1.252719
2	Services	0.4154754	0.2778976	0.1109746	0.1116213	0.0814478	1.185233
7	Services	0.5878499	0.1535348	0.0006138	0.1822722	0.0231219	1.123097
6	Services	0.3184271	0.2718329	0.2489164	0.0520256	0.1037725	1.128411
8	Retail Trade	0.3603968	0.1876330	0.0854220	0.0934442	0.0894848	1.119306

wide_topics[,c(2,6,9,10,11,12,13)] %>% filter(industry!="Finance") %>% arrange(desc(dist)) %>% mutate(id=1:n())%>% select(id,everything()) %>% slice(5) %>% html_df()

id	industry	Investment	Real Estate	Service2	Services	Utility	dist
5	Utilities	0.1768971	0.1143861	0.2481198	0.4017117	0.053171	1.144542

train <- wide_topics[!is.na(wide_topics$industry),]
label <- wide_topics[is.na(wide_topics$industry),]

library(caret)
tout <- readRDS('../../Data/corp_knn.rds')

tout

k-Nearest Neighbors 

5804 samples
  10 predictor
   9 classes: 'Agriculture', 'Construction', 'Finance', 'Manufacturing', 'Mining', 'Retail Trade', 'Services', 'Utilities', 'Wholesale Trade' 

No pre-processing
Resampling: Cross-Validated (10 fold) 
Summary of sample sizes: 5226, 5222, 5223, 5224, 5223, 5226, ... 
Resampling results across tuning parameters:

  k   Accuracy   Kappa    
   1  0.6922669  0.6037548
   2  0.6883222  0.5984635
   3  0.7219205  0.6397779
   4  0.7305403  0.6495724
   5  0.7374387  0.6581581
   6  0.7384702  0.6592123
   7  0.7460449  0.6686815
   8  0.7505306  0.6741651
   9  0.7515604  0.6753179
  10  0.7512102  0.6749574
  11  0.7489795  0.6718804
  12  0.7491537  0.6719035
  13  0.7525919  0.6764543
  14  0.7508766  0.6741010
  15  0.7529349  0.6766597
  16  0.7506983  0.6737148
  17  0.7500110  0.6727821
  18  0.7488041  0.6711643
  19  0.7494908  0.6718556
  20  0.7496676  0.6719961

Accuracy was used to select the optimal model using the largest value.
The final value used for the model was k = 15.

ggplot(tout$results, aes(x=k, y=Accuracy)) +
  geom_line() + 
  geom_ribbon(aes(ymin=Accuracy - AccuracySD*1.96,
                  ymax=Accuracy + AccuracySD*1.96), alpha=0.2) + 
  geom_vline(xintercept=15, color="blue") + 
  xlab("k, optimal = 15")

label$industry_pred <- predict(tout,
                              label)
label[,c("document",
         "industry_pred")] %>%
  head %>% html_df

document	industry_pred
0000817473-14-000010.txt	Finance
0000820027-14-000025.txt	Finance
0000837465-14-000002.txt	Manufacturing
0000837919-14-000002.txt	Finance
0000891092-14-000570.txt	Finance
0000891092-14-002078.txt	Finance

ts_wt <- wide_nodup %>% left_join(label[,c("document","industry_pred")])

ts_wt <- ts_wt %>%
  mutate(tsne1 = tsne_data$Y[, 1], tsne2 = tsne_data$Y[, 2])

# Force consistent factor values
inds <- unique(ts_wt$industry)
ts_wt$industry <- factor(ts_wt$industry, inds)
ts_wt$industry_pred <- factor(ts_wt$industry_pred, inds)

# Replicate default ggplot colors
ggplotColours <- function(n = 6, h = c(0, 360) + 15){
  if ((diff(h) %% 360) < 1) h[2] <- h[2] - 360/n
  hcl(h = (seq(h[1], h[2], length = n)), c = 100, l = 65)
}

ggplot() +
  scale_shape_identity() + # Allow for more plot point options
  geom_point(data=ts_wt[!is.na(ts_wt$industry),],
             aes(x=tsne1, y=tsne2, color=industry, shape=1), size=1) + 
  geom_point(data=ts_wt[!is.na(ts_wt$industry_pred),], aes(x=tsne1, y=tsne2,
             fill=industry_pred, shape=23, stroke=0.5), size=2) +
  guides(fill = "none") + 
  scale_color_manual(values=ggplotColours(n = 9), labels=inds, drop=FALSE) + 
  scale_fill_manual(values=ggplotColours(n = 9), labels=inds, drop=FALSE)

clusters <- kmeans(mat, 20)
clusters$cluster %>% head()

[1] 10  3 17  1  2 10

wide_nodup$kmean2 <- clusters$cluster[-dups]
ggplot(wide_nodup, aes(x = tsne1, y = tsne2, colour = factor(kmean2))) + 
    geom_point(alpha = 0.3) + theme_bw()

ggplot(wide_nodup, aes(x=kmean2)) + geom_bar() + facet_wrap(~factor(industry))

ggplot(wide_nodup, aes(x=tsne1, y=tsne2, color=factor(kmean2))) + geom_point() +
  facet_wrap(~factor(industry))