R 机器学习快速启动指南（三）
原文：
annas-archive.org/md5/13dee0bdb6445b090ed411f424dc82f4译者：飞龙
第六章：可视化欧洲联盟的经济问题
现在，我们将继续探讨第二个问题，即在不同国家，尤其是在欧盟内部检测宏观经济失衡的问题。在我看来，这失败了，并加剧了金融危机的损害。
在本章中，我们将探讨我们的问题，并为那些被确定为遭受宏观经济失衡的国家创建聚类。我们将涵盖以下主题：
-
欧洲联盟经济问题的概述
-
根据宏观经济失衡对国家进行聚类
各国经济问题的概述
2008 年的全球金融危机始于银行层面的危机，关闭了资金流动并增加了金融不稳定。在前几个月，主要影响的是金融机构和中介机构。然而，之后它达到了宏观经济层面，威胁到整个国家的稳定。
欧洲联盟、国际货币基金组织（IMF）和欧洲中央银行（ECB）采取了紧急措施以避免某些国家的破产。因此，希腊在 2010 年获得了 1100 亿欧元来偿还其公共债务。这笔贷款是进一步支出的开始。爱尔兰在 2017 年 11 月获得了 870 亿欧元，目的相同。随后，2011 年在葡萄牙发现了问题，希腊又收到了一笔额外的贷款。西班牙和塞浦路斯分别在 2012 年和 2013 年获得了资金。
当发生金融危机或出现新问题时，人们常说这是必然会发生的事情，或者每个人都预测到了。在金融危机之前，宏观经济失衡一直在稳步增加，但可能因为从 2002 年开始的全球经济增长而被隐藏。信用评级机构也无法检测到这些宏观经济失衡。
理解信用评级
信用评级可以描述为一个类别或等级，它代表了一个政府偿还其债务的能力和意愿，并按时全额偿还。与信用评级相关联的是一个潜在的违约概率。
信用评级对国家来说非常重要，因为它们决定了资本市场的融资成本，并使他们能够获得衍生品和贷款合同。
近年来，评级数量及其重要性显著增加。在一个更加全球化和复杂的世界中，投资者需要标准指标来比较发行人之间的信用质量，即使他们属于非常不同的国家。尽管信用评级机构标准普尔认为这并不是一门精确的科学，但信用评级根据各国的偿债能力进行排名。它们被认为是投资者和国际组织的重要参考。在下表中，描述了不同的等级及其含义：
信用评级机构的作用
依赖信用评级赋予了信用评级机构高水平的影響力和政治权力。尽管有大量的信用评级机构，但只有其中三家收集了超过 90%的市场份额。这些是穆迪、标准普尔和惠誉。
信用评级机构因在 1990 年代末的亚洲和俄罗斯金融危机以及 2008 年爆发的近期全球金融危机等重要危机中的表现而受到批评。在这些事件中，信用评级机构加剧了不平衡，反应或预测违约事件的时间过长，然后过度反应，导致严重降级。
例如，在最近的金融危机中，欧洲经济体的评级平均下降了三个等级。希腊和意大利是受影响最严重的国家之一。
信用评级流程
评级过程中考虑了许多因素，包括政治、增长、外债、金融部门、公私部门结构、社会发展以及贸易。例如，标准普尔在六点量表上对国家的信用度进行五项关键因素的评分。然而，不同因素的权重以及这些因素的结合方式是未知的。此外，这些因素的重要性随时间变化。定性判断通常在评级分配过程中起着非常重要的作用。
对评级质量的担忧在过去二十年里激励了世界各地的研究人员。了解评级的决定因素对于减少对信用评级机构的依赖以及使信用评级的分配更加客观至关重要。
信用评级的复制主要基于不同的计量经济学模型，试图找到一组宏观经济变量来更好地解释外部评级。然而，信用评级研究广泛忽视了定性和更为主观的信息。这些附加信息可以补充当前的定量宏观经济信息，并提高信用评级模型的性能。
根据宏观经济不平衡对国家进行聚类
在本节中，我们将开发一个无监督模型来直观地检测国家的宏观经济问题，并更深入地了解信用评级的主要驱动因素。我们将首先创建一个具有宏观经济问题的国家集群。在下一章中，我们将继续基于这些集群预测信用评级。
在本章中，我试图重复前几章中的代码。
让我们开始吧！
数据收集
与之前的模型一样，我们需要收集尽可能多的数据。首先，我们需要国家的宏观经济指标来分析宏观经济不平衡并预测主权评级。
下载和查看数据
我们将使用wbstats包，它提供了对从世界银行 API 获取的所有信息的结构化访问，包括所有年度、季度和月度数据。在这个包中，wb_cachelist函数提供了一个可用国家、指标和其他相关信息的快照，如下所示：
library(wbstats)
str(wb_cachelist, max.level = 1)
## List of 7
## $ countries :'data.frame': 304 obs. of 18 variables:
## $ indicators :'data.frame': 16978 obs. of 7 variables:
## $ sources :'data.frame': 43 obs. of 8 variables:
## $ datacatalog:'data.frame': 238 obs. of 29 variables:
## $ topics :'data.frame': 21 obs. of 3 variables:
## $ income :'data.frame': 7 obs. of 3 variables:
## $ lending :'data.frame': 4 obs. of 3 variables:
超过 16,000 个变量适用于 300 多个国家。默认语言是英语。
要搜索不同的指标，请使用wbsearch函数。例如，我们可以搜索包含国内生产总值（GDP）的所有指标，如下所示：
new_cache <- wbcache()
gdp_vars <- wbsearch(pattern = "gdp")
print(gdp_vars[1:5,])
## indicatorID
## 2 XGDP.56.FSGOV.FDINSTADM.FFD
## 3 XGDP.23.FSGOV.FDINSTADM.FFD
## 758 TG.VAL.TOTL.GG.ZS
## 759 TG.VAL.TOTL.GD.ZS
## 760 TG.VAL.TOTL.GD.PP.ZS
## indicator
## Government expenditure in tertiary institutions as % of GDP (%)
## Government expenditure in secondary institutions education as % of GDP (%)
## Trade in goods (% of goods GDP)
## Merchandise trade (% of GDP)
## Trade (% of GDP, PPP)
您也可以下载数据并指定指标、开始日期和结束日期，如下所示：
stock_return <- wb(indicator = "GFDD.OM.02", startdate = 2000, enddate = 2017)
head(stock_return)
## iso3c date value indicatorID indicator
## 110 ARG 2016 31.1342 GFDD.OM.02 Stock market return (%, year- on-year)
## 111 ARG 2015 36.6400 GFDD.OM.02 Stock market return (%, year- on-year)
## 112 ARG 2014 103.1500 GFDD.OM.02 Stock market return (%, year- on-year)
## 113 ARG 2013 60.0070 GFDD.OM.02 Stock market return (%, year- on-year)
## 114 ARG 2012 -19.8370 GFDD.OM.02 Stock market return (%, year- on-year)
## 115 ARG 2011 22.1574 GFDD.OM.02 Stock market return (%, year- on-year)
## iso2c country
## 110 AR Argentina
## 111 AR Argentina
## 112 AR Argentina
## 113 AR Argentina
## 114 AR Argentina
## 115 AR Argentina
我们现在已下载了所有国家的股市回报率。前几个值属于阿根廷。也可以只为特定国家获取信息。让我们获取 2015 年至 2017 年美国和西班牙的总人口，如下所示：
population_stock <- wb(country = c("ES","US"),indicator = c("SP.POP.TOTL","GFDD.OM.02"), startdate = 2015, enddate = 2017)
head(population_stock)
## iso3c date value indicatorID indicator iso2c country
## 1 ESP 2017 46572028 SP.POP.TOTL Population, total ES Spain
## 2 ESP 2016 46484062 SP.POP.TOTL Population, total ES Spain
## 3 ESP 2015 46444832 SP.POP.TOTL Population, total ES Spain
## 4 USA 2017 325719178 SP.POP.TOTL Population, total US United States
## 5 USA 2016 323405935 SP.POP.TOTL Population, total US United States
## 6 USA 2015 321039839 SP.POP.TOTL Population, total US United States
变量总是按行获取，但这并不是展示信息的一种常见方式。使用return_wide = TRUE选项，如果获取了多个变量，它们将以列的形式展示。让我们看一下以下示例：
population_stock <- wb(country = c("ES","US"),
indicator = c("SP.POP.TOTL","GFDD.OM.02"), startdate = 2015, enddate =2017,return_wide = TRUE)
head(population_stock)
## iso3c date iso2c country GFDD.OM.02 SP.POP.TOTL
## 1 ESP 2015 ES Spain 1.96478 46444832
## 2 ESP 2016 ES Spain -18.17990 46484062
## 3 ESP 2017 ES Spain NA 46572028
## 4 USA 2015 US United States 6.71498 321039839
## 5 USA 2016 US United States 1.70065 323405935
## 6 USA 2017 US United States NA 325719178
如果一个指标的最新可用日期未知，mrv参数代表最新值，并取一个整数，对应于所需的最新的值的数量。让我们显示美国股市回报率的 10 个最新值，如下所示：
wb(country = c("US"),indicator = "GFDD.OM.02", mrv=10,return_wide = TRUE)
## iso3c date iso2c country GFDD.OM.02
## 1 USA 2007 US United States 12.72230
## 2 USA 2008 US United States -17.40770
## 3 USA 2009 US United States -22.29390
## 4 USA 2010 US United States 20.24370
## 5 USA 2011 US United States 11.19970
## 6 USA 2012 US United States 8.81282
## 7 USA 2013 US United States 19.17170
## 8 USA 2014 US United States 17.49470
## 9 USA 2015 US United States 6.71498
## 10 USA 2016 US United States 1.70065
默认日期格式对排序或绘图没有用。POSIXct参数添加了date_ct和granularity列，使这些任务变得容易得多。
让我们以查看美国军事费用随时间演变的例子为例，如下所示：
library(ggplot2)
military_exp <- wb(country = c("US"),indicator = "MS.MIL.XPND.GD.ZS", POSIXct = TRUE)
ggplot() + theme_bw() +
geom_line(aes(y = value, x = date_ct), size=1.5, data = military_exp,
stat="identity") +
theme(legend.position="bottom", legend.direction="horizontal",
legend.title = element_blank()) +
labs(x="Year", y="Expenses as %GDP") +
ggtitle("US Military expenses %GDP")
您将得到以下输出图：
我们看到下载信息非常容易。然而，为了开发我们的模型，我们不会使用所有数据。现在，让我们定义我们需要多少变量、国家和历史数据。
数据简化
让我们使用以下代码查看世界银行数据库中的国家列表：
countries<-wb_cachelist$countries
尽管列表很长，但一些国家属于几个国家的联盟，例如阿拉伯世界或欧盟。我们需要对这些国家进行选择。
让我们查看以下可用的指标列表：
indicators<-wb_cachelist$indicators
指标列表甚至更大（有超过 16,000 个条目），但我们将选择其中最重要的，如下面的片段所示。这些指标是通过主要信用评级机构提供的评级方法指南获得的：
relevant_indicators<-c('NYGDPMKTPKDZ','FB.BNK.CAPA.ZS','GFDD.OI.01','GFDD.EI.07','GFDD.SI.04','GFDD.OI.02','GFDD.EI.02','FD.RES.LIQU.AS.ZS','FB.AST.NPER.ZS','GFDD.SI.05','GFDD.EI.05','GFDD.EI.09','GFDD.EI.06','GFDD.EI.10','GFDD.SI.01','FM.LBL.BMNY.GD.ZS','FM.LBL.BMNY.ZG','FS.AST.CGOV.GD.ZS','CC.EST','GFDD.EI.08','BN.CAB.XOKA.GD.ZS','IC.CRD.INFO.XQ','FS.AST.DOMS.GD.ZS','NE.EXP.GNFS.KD.ZG','NE.RSB.GNFS.ZS','GFDD.DI.08','NY.GDP.MKTP.KD.ZG','NY.GDP.PCAP.CD','NY.GDP.PCAP.KD.ZG','NE.CON.GOVT.ZS','NE.CON.GOVT.KD.ZG','GE.EST','NY.GDS.TOTL.ZS','NE.GDI.FTOT.ZS','NE.GDI.FTOT.KD.ZG','NE.CON.PRVT.KD.ZG','NE.CON.PRVT.PC.KD.ZG','NE.IMP.GNFS.KD.ZG','NV.IND.TOTL.ZS','NV.IND.TOTL.KD.ZG','FP.CPI.TOTL.ZG','FR.INR.LNDP','CM.MKT.LCAP.GD.ZS','PV.EST','SP.POP.GROW','GFDD.SI.07','REER','RQ.EST','RL.EST','NV.SRV.TETC.ZS','NV.SRV.TETC.KD.ZG','DT.DOD.DSTC.ZS','DT.DOD.DSTC.IR.ZS','GFDD.OM.02','IC.LGL.CRED.XQ','TOTRESV','SL.UEM.TOTL.ZS','SL.UEM.1524.ZS','VA.EST','SP.POP.TOTL')
indicators<-indicators[indicators$indicatorID %in% relevant_indicators,]
head(indicators[,c("indicatorID","indicator")])
获得的指标列表如下：
## indicatorID
## 284 VA.EST
## 1793 TOTRESV
## 1986 PV.EST
## 2569 RQ.EST
## 2607 RL.EST
## 2636 REER
## indicator
## 284 Voice and Accountability: Estimate
## 1793 Total Reserves
## 1986 Political Stability and Absence of Violence/Terrorism: Estimate
## 2569 Regulatory Quality: Estimate
## 2607 Rule of Law: Estimate
## 2636 Real Effective Exchange Rate
让我们从 2000 年到 2018 年下载这些指标的历史数据。使用以下代码，这将花费大约五分钟：
macroeconomic_data<-wb(indicator = relevant_indicators,startdate = 2000, enddate = 2018,return_wide = TRUE,POSIXct = TRUE)
让我们如下获取我们数据的结构：
str(macroeconomic_data)
完整的输出如下：
## 'data.frame': 5903 obs. of 66 variables:
## $ iso3c : chr NA NA NA NA ...
## $ date : chr "2018" "2017" "2016" "2015" ...
## $ iso2c : chr "EA" "EA" "EA" "EA" ...
## $ country : chr "Euro Area" "Euro Area" "Euro Area" "Euro Area" ...
## $ date_ct : Date, format: "2018-01-01" "2017-01-01" ...
## $ granularity : chr "annual" "annual" "annual" "annual" ...
## $ BN.CAB.XOKA.GD.ZS : num NA NA NA NA NA NA NA NA NA NA ...
## $ CC.EST : num NA NA NA NA NA NA NA NA NA NA ...
## $ CM.MKT.LCAP.GD.ZS : num NA NA NA NA NA NA NA NA NA NA ...
## $ DT.DOD.DSTC.IR.ZS : num NA NA NA NA NA NA NA NA NA NA ...
## $ DT.DOD.DSTC.ZS : num NA NA NA NA NA NA NA NA NA NA ...
## $ FB.AST.NPER.ZS : num NA NA NA NA NA NA NA NA NA NA ...
## $ FB.BNK.CAPA.ZS : num NA NA NA NA NA NA NA NA NA NA ...
## $ FD.RES.LIQU.AS.ZS : num NA NA NA NA NA NA NA NA NA NA ...
## $ FM.LBL.BMNY.GD.ZS : num NA NA NA NA NA NA NA NA NA NA ...
## $ FM.LBL.BMNY.ZG : num NA NA NA NA NA NA NA NA NA NA ...
## $ FP.CPI.TOTL.ZG : num NA NA NA NA NA NA NA NA NA NA ...
## $ FR.INR.LNDP : num NA NA NA NA NA NA NA NA NA NA ...
## $ FS.AST.CGOV.GD.ZS : num NA NA NA NA NA NA NA NA NA NA ...
## $ FS.AST.DOMS.GD.ZS : num NA NA NA NA NA NA NA NA NA NA ...
## $ GE.EST : num NA NA NA NA NA NA NA NA NA NA ...
## $ GFDD.DI.08 : num NA NA NA NA NA NA NA NA NA NA ...
## $ GFDD.EI.02 : num NA NA NA NA NA NA NA NA NA NA ...
## $ GFDD.EI.05 : num NA NA NA NA NA NA NA NA NA NA ...
## $ GFDD.EI.06 : num NA NA NA NA NA NA NA NA NA NA ...
## $ GFDD.EI.07 : num NA NA NA NA NA NA NA NA NA NA ...
## $ GFDD.EI.08 : num NA NA NA NA NA NA NA NA NA NA ...
## $ GFDD.EI.09 : num NA NA NA NA NA NA NA NA NA NA ...
## $ GFDD.EI.10 : num NA NA NA NA NA NA NA NA NA NA ...
## $ GFDD.OI.01 : num NA NA NA NA NA NA NA NA NA NA ...
## $ GFDD.OI.02 : num NA NA NA NA NA NA NA NA NA NA ...
## $ GFDD.OM.02 : num NA NA NA NA NA NA NA NA NA NA ...
## $ GFDD.SI.01 : num NA NA NA NA NA NA NA NA NA NA ...
## $ GFDD.SI.04 : num NA NA NA NA NA NA NA NA NA NA ...
## $ GFDD.SI.05 : num NA NA NA NA NA NA NA NA NA NA ...
## $ GFDD.SI.07 : num NA NA NA NA NA NA NA NA NA NA ...
## $ IC.CRD.INFO.XQ : num NA NA NA NA NA NA NA NA NA NA ...
## $ IC.LGL.CRED.XQ : num NA NA NA NA NA NA NA NA NA NA ...
## $ NE.CON.GOVT.KD.ZG : num NA NA NA NA NA NA NA NA NA NA ...
## $ NE.CON.GOVT.ZS : num NA NA NA NA NA NA NA NA NA NA ...
## $ NE.CON.PRVT.KD.ZG : num NA NA NA NA NA NA NA NA NA NA ...
## $ NE.CON.PRVT.PC.KD.ZG: num NA NA NA NA NA NA NA NA NA NA ...
## $ NE.EXP.GNFS.KD.ZG : num NA NA NA NA NA NA NA NA NA NA ...
## $ NE.GDI.FTOT.KD.ZG : num NA NA NA NA NA NA NA NA NA NA ...
## $ NE.GDI.FTOT.ZS : num NA NA NA NA NA NA NA NA NA NA ...
## $ NE.IMP.GNFS.KD.ZG : num NA NA NA NA NA NA NA NA NA NA ...
## $ NE.RSB.GNFS.ZS : num NA NA NA NA NA NA NA NA NA NA ...
## $ NV.IND.TOTL.KD.ZG : num NA NA NA NA NA NA NA NA NA NA ...
## $ NV.IND.TOTL.ZS : num NA NA NA NA NA NA NA NA NA NA ...
## $ NV.SRV.TETC.KD.ZG : num NA NA NA NA NA NA NA NA NA NA ...
## $ NV.SRV.TETC.ZS : num NA NA NA NA NA NA NA NA NA NA ...
## $ NY.GDP.MKTP.KD.ZG : num NA NA NA NA NA NA NA NA NA NA ...
## $ NY.GDP.PCAP.CD : num NA NA NA NA NA NA NA NA NA NA ...
## $ NY.GDP.PCAP.KD.ZG : num NA NA NA NA NA NA NA NA NA NA ...
## $ NY.GDS.TOTL.ZS : num NA NA NA NA NA NA NA NA NA NA ...
## $ NYGDPMKTPKDZ : num 2.1 2.39 1.8 2.07 2.21 ...
## $ PV.EST : num NA NA NA NA NA NA NA NA NA NA ...
## $ REER : num NA NA NA NA NA NA NA NA NA NA ...
## $ RL.EST : num NA NA NA NA NA NA NA NA NA NA ...
## $ RQ.EST : num NA NA NA NA NA NA NA NA NA NA ...
## $ SL.UEM.1524.ZS : num NA NA NA NA NA NA NA NA NA NA ...
## $ SL.UEM.TOTL.ZS : num NA NA NA NA NA NA NA NA NA NA ...
## $ SP.POP.GROW : num NA NA NA NA NA NA NA NA NA NA ...
## $ SP.POP.TOTL : num NA NA NA NA NA NA NA NA NA NA ...
## $ TOTRESV : num NA NA NA NA NA NA NA NA NA NA ...
## $ VA.EST : num NA NA NA NA NA NA NA NA NA NA ...
如您所见，大多数变量在前面几行显示缺失值。让我们获取一个国家在前面的列表中出现的次数列表：
head(table(macroeconomic_data$country))
##
## Advanced economies Advanced Economies Afghanistan
## 4 19 19
## Africa Albania Algeria
## 12 19 19
该数据包含的国家总数如下：
length(table(macroeconomic_data$country))
## [1] 330
我们可以如下探索我们的缺失值：
library(DataExplorer)
plot_missing(macroeconomic_data)
我们将获得以下输出：
获得的图表不太美观，因为变量数量众多。虽然信息容易获取，但许多变量显示大量缺失值。这是因为信息未更新，且最近几年的数据不可用。
要查看结果，您应手动选择不同官方来源的相关变量并在 Excel 文件中合并它们，该文件将用于我们的模型开发。首先，按照如下操作删除工作区中的所有对象：
rm(list=ls())
现在，将新的数据样本加载到 R 中，如下：
macroeconomic_data <- read.csv("SovereignData.csv",sep=";",header=TRUE,stringsAsFactors = FALSE,dec=",")
该文件包含以下信息：
str(macroeconomic_data)
## 'data.frame': 224 obs. of 14 variables:
## $ Year : int 2010 2011 2012 2013 2014 2015 2016 2017 2010 2011 ...
## $ CountryISO : chr "AT" "AT" "AT" "AT" ...
## $ CountryName : chr "Austria" "Austria" "Austria" "Austria" ...
## $ CARA : num 2.92 1.57 1.5 1.95 2.39 ...
## $ DCPI : num 1.69 3.54 2.57 2.12 1.47 ...
## $ DGDP : num 1.828 2.935 0.662 0.008 0.923 ...
## $ ILMA : num 22284 25161 27211 23290 24941 ...
## $ PSBR : num -4.44 -2.55 -2.19 -1.95 -2.68 ...
## $ PUDP : num 82.4 82.2 81.7 81 83.7 ...
## $ TDRA : num -0.469 -1.168 -0.992 -0.313 0.232 ...
## $ YPCA : num 0.0467 0.0511 0.0483 0.0506 0.0519 ...
## $ MEANWGI : num 1.55 1.48 1.52 1.54 1.52 ...
## $ YearsFromLastDefault: int 21 22 23 24 25 26 27 28 21 22 ...
## $ RatingMayT1 : int 6 5 5 5 5 5 5 5 5 5 ...
如str函数所示，数据样本包含224个观测值和以下14个变量：
-
Year: 年末的参考日期 -
CountryISO: 国际标准国家代码 -
CountryName: 国家的名称 -
CARA: 一个国家的经常账户余额 -
DCPI: 消费者价格增长率 -
DGDP: GDP 增长率 -
ILMA: 国际储备 -
PSBR: 预算平衡占 GDP 的百分比 -
PUDP: 公共债务占 GDP 的百分比 -
TDRA: 总进口减去总出口的商品和服务占 GDP 的百分比 -
YPCA: 人均 GDP -
MEANWGI: 六个全球治理指标的平均值 -
YearsFromLastDefault: 自上次国家违约以来经过的时间 -
RatingMayT1: 宏观经济信息一年后国家的信用评级
您可以从以下来源获取更多数据：
所考虑的所有变量在先前的研究或由信用评级机构（CRAs）使用中已被广泛采用。在我们的变量列表中，我们使用了世界银行的多种治理指标，包括MEANWGI变量。该变量是通过世界银行提供的六个不同指标计算得出的简单平均值。这些指标衡量的话题包括政治稳定性、暴力水平和对腐败的控制。世界银行根据不同非政府组织、国际组织或私营部门公司的调查来计算这些综合指标。每个指标的价值越高，一个国家在该指标中的质量或评分就越高。
本章剩余部分将使用的数据集收集了多年来每个国家在五月底的信用评级。宏观经济信息代表了每个国家上一年年末的经济状况。各国在年初的第一季度或第二季度发布这些信息，因此每年五月份获取的最新宏观经济信息属于上一年。
研究数据
我们的首要目标将是根据各国的宏观经济状况检测国家集群，并预测这些国家的评级。然后，我们将评估欧洲委员会发布的国家报告中的定性信息是否有助于预测国家评级。这些报告通常在每年的前几个月发布，因此它们应该在五月的信用评级中也有影响力。
再次，我们必须从分析以下数据集开始：
head(macroeconomic_data)
## Year CountryISO CountryName CARA DCPI DGDP ILMA PSBR PUDP TDRA
## 1 2010 AT Austria 2.923 1.694 1.828 22284 -4.440 82.401 -0.469
## 2 2011 AT Austria 1.575 3.542 2.935 25161 -2.554 82.193 -1.168
## 3 2012 AT Austria 1.499 2.569 0.662 27211 -2.189 81.667 -0.992
## 4 2013 AT Austria 1.947 2.118 0.008 23290 -1.950 81.016 -0.313
## 5 2014 AT Austria 2.389 1.468 0.923 24941 -2.683 83.714 0.232
## 6 2015 AT Austria 1.932 0.810 1.073 22236 -1.033 84.297 0.609
## YPCA MEANWGI YearsFromLastDefault RatingMayT1
## 1 0.0466532 1.552740 21 6
## 2 0.0510888 1.477206 22 5
## 3 0.0483316 1.515706 23 5
## 4 0.0505913 1.540124 24 5
## 5 0.0519292 1.521760 25 5
## 6 0.0447088 1.471732 26 5
该数据集中没有缺失值。让我们按以下方式探索我们的数据：
library(funModeling)
library(DataExplorer)
info<-df_status(macroeconomic_data)
##variable q_zeros p_zeros q_na p_na q_inf p_inf type
## Year 0 0 0 0 0 0 integer
## CountryISO 0 0 0 0 0 0 character
## CountryName 0 0 0 0 0 0 character
## CARA 0 0 0 0 0 0 numeric
## DCPI 0 0 0 0 0 0 numeric
## DGDP 0 0 0 0 0 0 numeric
## ILMA 0 0 0 0 0 0 numeric
## PSBR 0 0 0 0 0 0 numeric
## PUDP 0 0 0 0 0 0 numeric
## TDRA 0 0 0 0 0 0 numeric
## YPCA 0 0 0 0 0 0 numeric
## MEANWGI 0 0 0 0 0 0 numeric
## YearsFromLastDefault 0 0 0 0 0 0 integer
## RatingMayT1 0 0 0 0 0 0 integer
不同国家的宏观经济信息从 2010 年到 2017 年可用，如下所示：
table(macroeconomic_data$Year)
## 2010 2011 2012 2013 2014 2015 2016 2017
## 28 28 28 28 28 28 28 28
以下是需要考虑的 28 个国家的列表：
unique(macroeconomic_data$CountryName)
## [1] "Austria" "Belgium" "Bulgaria" "Croatia"
## [5] "Cyprus" "Czech Republic" "Denmark" "Estonia"
## [9] "Finland" "France" "Germany" "Greece"
## [13] "Hungary" "Ireland" "Italy" "Latvia"
## [17] "Lithuania" "Luxembourg" "Malta" "Netherlands"
## [21] "Poland" "Portugal" "Romania" "Slovakia"
## [25] "Slovenia" "Spain" "Sweden" "United Kingdom"
获取目标变量
目标变量将是每个国家的信用评级。让我们用以下代码查看目标变量的不同值：
unique(macroeconomic_data$RatingMayT1)
## [1] 6 5 3 2 4 1
如我们所见，该量表只包含六个不同的等级，它们是数字的，而不是字母的。以下代码创建了一个表格来显示每个评级类别的分配数字：
RatingSP <- c('AAA','AA+','AA','AA-','A+','A','A-','BBB+','BBB','BBB-','BB+','BB','BB-','B+','B','B-','CCC+','CCC','CCC-','CC','C','D','DD','SD')
Rating_num <- c('6','5','5','5','4','4','4','3','3','3','2','2','2','1','1','1','1','1','1','1','1','1','1','1')
mapping<-data.frame(RatingSP, Rating_num)
rm(RatingSP,Rating_num)
print(mapping)
以下映射表将被用于将所有可能的信用评级水平减少到仅六个类别：
## RatingSP Rating_num
## 1 AAA 6
## 2 AA+ 5
## 3 AA 5
## 4 AA- 5
## 5 A+ 4
## 6 A 4
## 7 A- 4
## 8 BBB+ 3
## 9 BBB 3
## 10 BBB- 3
## 11 BB+ 2
## 12 BB 2
## 13 BB- 2
## 14 B+ 1
## 15 B 1
## 16 B- 1
## 17 CCC+ 1
## 18 CCC 1
## 19 CCC- 1
## 20 CC 1
## 21 C 1
## 22 D 1
## 23 DD 1
## 24 SD 1
因此，6的值与最高评级等级（AAA）相关联，而1则对应最低评级等级。这种映射是为了减少不同标签的数量，并获得一个更细粒度的目标变量。
让我们用以下代码查看目标变量的分布情况：
tab<-table(macroeconomic_data$RatingMayT1)
barplot(tab,xlab="Rating",ylab="Count",border="blue",col="blue")
你将看到一个如下所示的图表：
获取信用质量
2018 年 5 月，每个欧洲国家的信用评级如下：
with(macroeconomic_data[macroeconomic_data$Year==2017,], table(CountryName,RatingMayT1))
每个国家的信用质量表如下：
## RatingMayT1
## CountryName 1 2 3 4 5 6
## Austria 0 0 0 0 1 0
## Belgium 0 0 0 0 1 0
## Bulgaria 0 1 0 0 0 0
## Croatia 0 1 0 0 0 0
## Cyprus 0 1 0 0 0 0
## Czech Republic 0 0 0 0 1 0
## Denmark 0 0 0 0 0 1
## Estonia 0 0 0 0 1 0
## Finland 0 0 0 0 1 0
## France 0 0 0 0 1 0
## Germany 0 0 0 0 0 1
## Greece 1 0 0 0 0 0
## Hungary 0 0 1 0 0 0
## Ireland 0 0 0 1 0 0
## Italy 0 0 1 0 0 0
## Latvia 0 0 0 1 0 0
## Lithuania 0 0 0 1 0 0
## Luxembourg 0 0 0 0 0 1
## Malta 0 0 0 1 0 0
## Netherlands 0 0 0 0 0 1
## Poland 0 0 1 0 0 0
## Portugal 0 1 0 0 0 0
## Romania 0 0 1 0 0 0
## Slovakia 0 0 0 1 0 0
## Slovenia 0 0 0 1 0 0
## Spain 0 0 1 0 0 0
## Sweden 0 0 0 0 0 1
## United Kingdom 0 0 0 0 1 0
丹麦、德国、卢森堡、荷兰和瑞典显示出最高的信用质量。另一方面，根据标准普尔的数据，希腊是欧盟中偿债能力最低的国家。
在地图上显示信用评级
使用rworldmap包，我们可以创建一个地图并按国家显示信用评级。让我们关注最后可用的信用评级（2018 年 5 月），如下所示：
macro2017<-macroeconomic_data[macroeconomic_data$Year==2017,]
library(rworldmap)
Map <- joinCountryData2Map(macro2017, joinCode = "ISO2",
nameJoinColumn = "CountryISO")
可以快速总结如下：
## 28 codes from your data successfully matched countries in the map
## 0 codes from your data failed to match with a country code in the map
## 215 codes from the map weren't represented in your data
指定轴的极限，如下所示：
mapCountryData(Map, catMethod = "categorical", missingCountryCol = gray(.8), xlim = c(-20, 59),ylim = c(35, 71),asp = 1)
将会显示如下非常漂亮的输出：
执行数据的描述性分析
我们现在应该进行描述性分析，就像我们在前面的章节中所做的那样，使用以下代码：
library(fBasics)library(DataExplorer)
descriptives_num<-as.data.frame(t(basicStats(macroeconomic_data[,4:13])))
head(descriptives_num)
获得的输出如下：
## nobs NAs Minimum Maximum 1\. Quartile 3\. Quartile Mean
## CARA 224 0 -11.354 12.638 -1.44275 4.15000 1.004080
## DCPI 224 0 -2.097 6.113 0.27950 2.44350 1.418147
## DGDP 224 0 -9.168 25.492 0.87175 3.17450 1.946969
## ILMA 224 0 207.470 248856.000 3099.44000 59386.00000 45290.293839
## PSBR 224 0 -32.055 3.931 -4.32325 -0.94425 -3.110054
## PUDP 224 0 6.067 181.147 40.68500 86.39600 68.819321
## Median Sum SE Mean LCL Mean UCL Mean
## CARA 0.6140 224.914 0.284204 0.444012 1.564149
## DCPI 1.3000 317.665 0.103081 1.215009 1.621285
## DGDP 1.8995 436.121 0.195218 1.562261 2.331676
## ILMA 23359.0000 10145025.820 3712.287504 37974.641205 52605.946474
## PSBR -2.5980 -696.652 0.239382 -3.581794 -2.638313
## PUDP 65.0070 15415.528 2.417700 64.054860 73.583783
## Variance Stdev Skewness Kurtosis
## CARA 1.809288e+01 4.253572 0.057613 0.070869
## DCPI 2.380167e+00 1.542779 0.324855 -0.127702
## DGDP 8.536634e+00 2.921752 1.873828 18.872568
## ILMA 3.086962e+09 55560.431839 1.508513 1.483938
## PSBR 1.283607e+01 3.582746 -2.870375 18.134018
## PUDP 1.309341e+03 36.184813 0.792637 0.645762
按如下方式绘制这些变量：
plot_histogram(macroeconomic_data[,4:13])
下面的图表如下所示：
一旦我们研究过我们的数据集，我们就可以开始处理我们的数据了。保存工作区的一个备份。如果出了问题，总是可以加载备份并继续使用代码，如下所示：
rm(list=setdiff(ls(), c("macroeconomic_data")))
save.image("Backup1.RData")
检测宏观经济失衡
在本节中，我们将使用称为自组织映射（SOM）的技术，根据宏观经济变量对欧洲国家进行分类。
自组织映射技术
识别国家之间的相似性和差异性对于检测失衡并采取纠正措施以避免金融危机的传播至关重要。SOM 是一种非常适用于探索性数据分析（EDA）的无监督神经网络。SOM 的无监督性质是由于网络能够学习数据中的模式，而不使用目标变量。因此，网络中包含的不同神经元通过输入数据自我组织。以下图是 SOM 最常见的表示形式：
如前图所示，存在两个主要的神经元层：一个输入层和一个输出层。输入层是我们发送到网络的输入数据，即要分析的数据。输出层是一个二维地图，其中将放置所有我们的数据。输入层中的数据通过突触权重连接到输出层的所有神经元。这样，输入层每个神经元提供的信息就传递到了输出层的所有神经元。
网络试图找到输出层中与输入层神经元值最相似的突触权重神经元。最相似的权重在输入数据和输出层中创建的原型权重之间具有最低的距离。这些权重最初是随机创建的。一旦计算出距离，就会识别出获胜神经元。这个神经元是权重与输入集之间差异最小或欧几里得距离最短的神经元。
一旦识别出获胜神经元，其权重将使用学习规则进行调整，以将这些权重代理到使神经元获胜的输入模式。这样，权重最接近输入模式的神经元被更新，以变得更加接近。因此，输出层中的每个神经元都专门化于这个输入模式。
获胜神经元不是唯一更新其权重的神经元。还使用邻域函数来更新获胜神经元及其邻域神经元的权重，以定位相似的模式。随着模型迭代次数的增加，这个邻域半径会减小，以实现每个神经元的更好专业化。
训练 SOM
现在是时候使用 R 来训练 SOM 了。为此，我们将使用kohonen包。这个包的名字来源于 Teuvo Kohonen，他首先介绍了这个算法。
由于 SOM 主要基于欧几里得距离，建议对变量进行标准化。我们将使用caret包来完成这项工作，如下所示：
library(kohonen)
library(caret)
preprocess <- preProcess(macroeconomic_data[,4:13], method=c("center", "scale"))
print(preprocess)
## Created from 224 samples and 10 variables
##
## Pre-processing:
## - centered (10)
## - ignored (0)
## - scaled (10)
这将创建一个包含转换变量的数据框。然后，我们将国家名称和评分添加到这个转换数据框中，如下所示：
macroeconomic_data_trn <- cbind(macroeconomic_data[,c(1:3,14)],predict(preprocess, macroeconomic_data[,4:13]))
现在，是时候训练映射了。网络的第一层将有十个输入模式或十个宏观经济变量。我们的输出层将被固定为一个 6×5 大小的二维映射。这个大小是在多次测试后获得的，如下所示：
set.seed(1234)
som_grid <- somgrid(xdim = 6, ydim=5, topo="hexagonal")
som_model <- som(as.matrix(macroeconomic_data_trn[,5:14]),grid=som_grid,rlen=800,alpha=c(0.1,0.01), keep.data = TRUE )
经过 800 次迭代后，训练过程结束。可以从训练好的映射中获取不同的图表。
让我们使用以下代码显示训练过程中到最近码本向量的平均距离：
plot(som_model, type = "changes")
这将绘制以下图表：
我们可以看到，训练过程中的误差，即测量到最近单元的平均距离，在迭代过程中是逐渐减少的。如果迭代次数足够高，算法就会收敛。
此外，还可以可视化映射并计算每个神经元中分类的国家数量：
plot(som_model, type = "counts", main="Node Counts")
节点计数的图形表示可以在以下图表中看到：
在图中，根据分类的国家数量，用不同颜色标记了 30 个单元格（大小为 6x5）。灰色单元格表示在这个单元格中没有国家被分类。单元格数量越多，每个神经元中的国家数量就越少，从而获得更细粒度的地图。
在下面的图中，显示了映射到地图中每个单元的国家之间的平均距离。距离越小，表示国家由代码簿向量表示得越好，如下所示：
plot(som_model, type = "quality", main="Node Quality/Distance")
节点质量图如下：
地图上的权重向量也可以被可视化。这有助于揭示不同国家分布中的模式。
节点权重向量或代码由用于生成 SOM 的原变量的归一化值组成：
plot(som_model, type = "codes")
前面代码的输出显示在以下图表中：
从前面的图表中，可以提取出被放置在每个单元格中的国家的特征。还可以可视化单个变量在地图上的分布情况。例如，让我们看一下MEANWGI变量的分布，如下所示：
plot(som_model, type = "property", property = getCodes(som_model)[,'MEANWGI'], main="WorldBank Governance Indicators")
训练好的地图被绘制出来，并且根据这些地图上国家所在位置的变量值来着色单元格，如下所示：
根据前面的地图，该变量值较高的国家将位于地图的右上角。
现在，让我们使用以下代码，利用 2017 年 12 月的宏观经济信息来查看这些国家的位置：
Labels<-macroeconomic_data_trn[,c("CountryName","Year")]
Labels$CountryName<-ifelse(Labels$Year!=2017,"",as.character(Labels$CountryName))
plot(som_model, type = "mapping",label=Labels$CountryName)
将根据经济变量将国家放置在地图上的地图将显示如下：
在宏观经济变量中值相似的国家更靠近。因此，希腊、葡萄牙和塞浦路斯在地图上非常接近。这些国家在过去几年中由于金融危机而遇到了重要问题。德国、法国和英国等国家根据宏观经济信息也相似。最后，我们可以通过考虑训练地图的代码簿向量来在这个地图上创建不同的聚类或组。可以使用层次聚类来选择组数。以下代码将获得五个不同的国家组：
clusters <- cutree(hclust(dist(getCodes(som_model))), 5)
通过运行以下代码，可以在地图上可视化最近创建的组：
plot(som_model, type="codes", bgcol = clusters, main = "Clusters")
add.cluster.boundaries(som_model, clusters)
下面的地图被绘制出来：
前面的图表给出了以下信息：
-
组 1：葡萄牙、希腊和塞浦路斯等国家。
-
第二组：包括波兰、马耳他、保加利亚和西班牙等国家。总的来说，这是一个由西欧和其他经济状况较弱的国家组成的集群。这个组比之前的组状况更好。
-
第三组：欧盟中更富裕的国家，如英国、德国和法国。
-
第四组和第五组：这些组各有一个单元格，其中一个包含爱尔兰。爱尔兰几年前经历了严重的经济危机，但现在的情况非常不同。这种差异如此之大，以至于它不能被包含在任何其他国家的组中。
如往常一样，您现在可以使用以下代码备份您的 工作空间：
save.image("Backup2.RData")
使用这种简单的聚类方法，可以可视化欧盟每个国家的不同情况。
摘要
在本章中，我们介绍了不同欧洲国家经历的经济危机。我们获得了数据并对其进行了分析。然后，我们开发了一个可视化工具，同时使用不同的变量来比较各国。
在下一章中，我们将尝试创建一个评级模型。我们将根据各国的经济状况分配不同的分数。让我们看看我们是否能够预测不同国家的评级！
第七章：主权危机 - 自然语言处理和主题建模
继续探讨欧洲国家的经济问题检测，在本章中，我们将尝试使用定量和定性信息复制标准普尔提供的国家评级。
本章是一个有趣的实际案例应用，因为我们将使用一些基本的文本挖掘技术来复制标准普尔的信用评级。为此，我们将使用欧洲委员会为欧洲成员国发布的国家报告。
我们将执行文本挖掘过程，提取单词组合或有用的术语来预测主权评级。
本章将涵盖以下主题：
-
使用宏观经济信息预测国家评级
-
实现决策树
-
使用欧洲国家报告预测主权评级
使用宏观经济信息预测国家评级
在我们的聚类模型中，如第六章所述，使用自组织图可视化欧盟的经济问题，我们使用了所有可用数据。现在，为了训练一个能够预测主权评级的模型，我们需要将数据分为两个样本：训练样本和测试样本。
这对我们来说并不新鲜。当我们试图开发不同的模型来预测银行的失败时，我们使用了caTools包来分割数据，同时考虑我们的目标变量。
这里再次使用相同的程序：
library(caTools)
index = sample.split(macroeconomic_data$RatingMayT1, SplitRatio = .75)
train_macro<-subset(macroeconomic_data, index == TRUE)
test_macro<-subset(macroeconomic_data, index == FALSE)
现在，你可以打印以下语句：
print(paste("The number of observations in the train sample is: ",nrow(train_macro),sep=""))
## [1] "The number of observations in the train sample is: 168"
print(paste("The number of observations in the test sample is: ",nrow(test_macro),sep=""))
## [1] "The number of observations in the test sample is: 56"
因此，测试样本代表了总数据的 25%。此外，不同信用评级的相对比例在训练样本和测试样本中都得到了保留。
这两个样本都将进行标准化。再次使用caret包：
library(caret)
变换应仅考虑训练样本，然后应用于测试样本：
preprocess <- preProcess(train_macro[,4:13], method=c("center", "scale"))
print(preprocess)
## Created from 168 samples and 10 variables
##
## Pre-processing:
## - centered (10)
## - ignored (0)
## - scaled (10)
这里是两个额外的训练和测试样本，将原始变量替换为转换后的变量：
train_macro_trn <- cbind(train_macro[,c(1:3,14)],predict(preprocess, train_macro[,4:13]))
test_macro_trn <- cbind(test_macro[,c(1:3,14)],predict(preprocess, test_macro[,4:13]))
让我们看看变量是如何与目标变量（评级）相关的。为此，我们首先将目标变量转换为因子变量，然后我们将根据每个变量的类别创建不同的箱线图：
library(ggplot2)
variables<-colnames(train_macro_trn[,5:14])
train_macro_trn$RatingMayT1<-as.factor(train_macro_trn$RatingMayT1)
for (i in 5:14)
{
library(ggplot2)
theme_set(theme_classic())
var<-colnames(train_macro_trn)[i]
data_aux<-train_macro_trn[,c(var,"RatingMayT1")]
colnames(data_aux)<-c("variable","RatingMayT1")
g <- ggplot(data_aux, aes(RatingMayT1,variable))
plot(g + geom_boxplot(varwidth=T, fill="plum") +
labs(title="Box plot",
subtitle=var,
x="Rating Number",
y=var))
}
以下是以一个国家的经常账户余额（CARA）变量为输出的结果：
以下是以消费者价格增长率（DCPI）变量为输出的结果：
以下是以 GDP 增长（DGPD）变量为输出的结果：
以下是以国际储备（ILMA）变量为输出的结果：
以下是以六个全球治理指标平均值（MEANWGI）变量为输出的结果：
以下是以 GDP（国内生产总值）百分比表示的预算平衡（PSBR）变量的输出：
以下是对公共债务比率（PUDP）变量的输出：
以下是对商品和服务总进口减去总出口占 GDP 百分比（TDRA）变量的输出：
以下是对 GDP 人均值（YPCA）变量的输出：
这些图表有助于观察数据中的某些模式以及变量如何帮助预测信用评级。
例如，对于公共债务（GDP 百分比）的箱形图，或PUDP，显示评级最低的国家在这个变量上的平均值较高。
在以下代码中，让我们使用上一个图表，但这次提供更多关于每个评级类别中包含的国家细节：
library(dplyr)
means_YPCA <- train_macro_trn %>% group_by(RatingMayT1) %>%
summarise(YPCA = mean(YPCA))
ggplot(train_macro_trn, aes(x = RatingMayT1, y = YPCA, color = RatingMayT1, fill = RatingMayT1)) +
geom_bar(data = means_YPCA, stat = "identity", alpha = .3) + ggrepel::geom_text_repel(aes(label = CountryName), color = "black", size = 2.5, segment.color = "grey") + geom_point() + guides(color = "none", fill = "none") + theme_bw() + labs( x = "Rating", y = "GDP per capita")
以下是对 YPCA 的更详细图表：
让我们看看另一个替代方案：
library(ggplot2)
theme_set(theme_classic())
ggplot(train_macro_trn, aes((MEANWGI))) + geom_density(aes(fill=factor(RatingMayT1)),alpha=0.8) +
labs(title="Density plot",
subtitle="Mean of the Worldwide Governance Indicators",
x=paste("MeanWGI",sep=''),
fill="RatingNum")
这是 MeanWGI 的密度图：
让我们为CARA实现一个Density plot：
ggplot(train_macro_trn, aes((CARA))) + geom_density(aes(fill=factor(RatingMayT1)),alpha=0.8) + labs(title="Density plot",
subtitle="Current account balance/GDP",
x=paste("CARA",sep=''),
fill="RatingNum")
CARA 上的密度图输出将如下所示：
作为目标值，我们有一个取六个不同值的变量。在这个问题中，与上一个问题不同，我们预测的是失败的银行，无法计算六个。
为了评估每个变量预测信用评级的能力，我们可以计算每个变量与目标变量的相关性：
colnames(train_macro_trn)
## [1] "Year" "CountryISO" "CountryName"
## [4] "RatingMayT1" "CARA" "DCPI"
## [7] "DGDP" "ILMA" "PSBR"
## [10] "PUDP" "TDRA" "YPCA"
## [13] "MEANWGI" "YearsFromLastDefault"
variables<-colnames(train_macro_trn[,c(4:14)])
使用以下代码，我们首先创建我们训练样本的副本，然后将目标变量转换为数值格式。这是因为无法用非数值变量计算相关性。
然后我们将使用cor函数计算相关性：
aux<-train_macro_trn
aux$RatingMayT1<-as.numeric(as.character(aux$RatingMayT1))
# Correlation matrix
correlations<-cor(aux[, variables], use="pairwise", method="pearson")
correlations_with_Rating<-as.matrix(correlations[1,])
接下来，打印这些相关性：
print(correlations_with_Rating)
## [,1]
## RatingMayT1 1.0000000
## CARA 0.3938594
## DCPI 0.1517755
## DGDP 0.1167254
## ILMA 0.3130267
## PSBR 0.2783237
## PUDP -0.4172153
## TDRA 0.3854816
## YPCA 0.6491449
## MEANWGI 0.8024756
## YearsFromLastDefault 0.5132374
与信用评级最相关的变量是治理指标的平均值（MEANWGI），其次是人均 GDP（YPCA）。在这两种情况下，变量的值越高，偿付能力或信用评级值就越高。
另一方面，相关性最弱的变量是消费者价格变化（DCPI）。所有变量都有正相关，除了PUDP。这意味着一个国家的债务越高，信用评级就越低。
根据文献和信用评级机构提供的方法论指南，所有变量都与信用评级有预期的符号。
在这一点上，我们应该保存我们的工作空间并删除任何不必要的对象：
rm(list=setdiff(ls(), c("macroeconomic_data","train_macro","test_macro","correlations_with_Rating","train_macro_trn","test_macro_trn")))
save.image("Backup3.RData")
如所示，观察值和变量的数量与我们在第二章中获得的银行数据集大不相同，预测银行失败 - 数据收集。
让我们尝试一些算法来预测信用评级。具体来说，在下一节中，我们将训练一个决策树和一个有序逻辑回归。
实现决策树
在之前第五章的“测试随机森林模型”部分（预测银行失败 - 多变量分析）中，我们查看随机森林时，简要介绍了决策树。在决策树中，训练样本根据最显著的独立变量分成两个或更多同质集合。在决策树中，找到将数据分成不同类别的最佳变量。信息增益和基尼指数是找到这个变量的最常见方法。然后，数据递归分割，扩展树的叶节点，直到达到停止标准。
让我们看看如何在 R 中实现决策树以及这个算法如何预测信用评级。
决策树在rpart包中实现。此外，rpart.plot包将有助于稍后可视化我们的训练模型。通过以下步骤实现这些包：
library(rpart)
library(rpart.plot)
要创建一个树，我们将使用rpart函数。必须指定四个参数：
-
公式：在格式目标：
~ predictor1+predictor2+…+predictorN -
数据：指定数据框
-
方法：
class用于分类树或anova用于回归树 -
控制：控制树增长的可选参数
在我们的情况下，以下控制参数被指定：
-
maxdepth：指最终树中任何节点的最大深度。它定义了分割的数量，换句话说，树可以增长多少，考虑到根节点深度为 0。 -
复杂度参数（或
cp）：此参数也用于控制决策树的大小。此参数可以被认为是增加决策树增长或复杂性的最小增益。如果添加新节点到树中不增加我们的拟合度，算法将停止增长。
让我们训练模型。首先创建一个包含我们变量的列表：
variables<-names(train_macro_trn[,4:14])
print(variables)
## [1] "RatingMayT1" "CARA" "DCPI"
## [4] "DGDP" "ILMA" "PSBR"
## [7] "PUDP" "TDRA" "YPCA"
## [10] "MEANWGI" "YearsFromLastDefault"
现在，训练了一个决策树：
set.seed(1234)
DT<-rpart(formula = RatingMayT1 ~ ., data = train_macro_trn[,c(variables)], control=rpart.control(maxdepth=5,cp=0.001))
模型训练完成后，可以使用summary函数打印模型给出的所有信息，尽管这次由于输出量较大而没有打印：
#summary(DT)
现在让我们使用决策树预测信用评级，包括训练样本和测试样本：
DT_pr_train <- data.frame(cbind(train_macro_trn$CountryName,train_macro_trn$Year,train_macro_trn$RatingMayT1,predict(DT, newdata=train_macro_trn, type="class")))
colnames(DT_pr_train)<-c("Country","year","Observed","Predicted")
DT_pr_test <- data.frame(cbind(test_macro_trn$CountryName,test_macro_trn$Year,test_macro_trn$RatingMayT1,predict(DT, newdata=test_macro_trn, type="class")))
colnames(DT_pr_test)<-c("Country","year","Observed","Predicted")
这是训练样本的混淆表：
table(DT_pr_train$Observed,DT_pr_train$Predicted)
## 1 2 3 4 5 6
## 1 6 2 0 0 0 0
## 2 0 16 5 1 1 0
## 3 1 4 22 4 2 0
## 4 0 0 7 25 0 0
## 5 0 0 7 1 25 1
## 6 0 0 0 2 1 35
训练好的模型能够预测训练样本中几乎所有的信用评级。现在让我们打印其在测试样本中的准确率：
table(DT_pr_test$Observed,DT_pr_test$Predicted)
## 1 2 3 4 5 6
## 1 2 0 0 1 0 0
## 2 0 3 5 0 0 0
## 3 0 1 8 2 0 0
## 4 0 0 1 8 1 0
## 5 0 0 2 2 7 1
## 6 0 0 0 1 1 10
为了评估决策树的准确性，我们可以计算不同的指标。具体来说，我们可以计算实际评级值与预测值之间的差异。通过计算这些差异，我们可以测量我们的模型在哪个等级上与实际评级水平不同。
为此，创建了一个函数：
model_assessment<-function(data,model)
{
data$Observed<-as.numeric(as.character(data$Observed))
data$Predicted<-as.numeric(as.character(data$Predicted))
data$df<-abs(as.numeric(data$Predicted)-as.numeric(data$Observed))
comparison<-as.data.frame(table(data$df))
comparison$perc<-comparison$Freq/nrow(data)
colnames(comparison)<-c("notche","N",paste("perc_",model,sep=''))
comparison$N<-NULL
comparison$cumulative<-cumsum(comparison[,ncol(comparison)])
return(comparison)
}
这里是不同的结果：
model_assessment(DT_pr_train,"DT")
## notche perc_DT cumulative
## 1 0 0.767857143 0.7678571
## 2 1 0.148809524 0.9166667
## 3 2 0.077380952 0.9940476
## 4 3 0.005952381 1.0000000
根据前面的表格，几乎 77%的国家被正确分类。另一方面，14.88%的国家的预测没有被正确分类，但与实际观察到的评级的差异只有一个等级。另一方面，7.74%的国家有一个错误的预测评级，并且这个预测与实际值相差两个等级，等等。
现在将相同的函数应用于测试样本：
model_assessment(DT_pr_test,"DT")
## notche perc_DT cumulative
## 1 0 0.67857143 0.6785714
## 2 1 0.25000000 0.9285714
## 3 2 0.05357143 0.9821429
## 4 3 0.01785714 1.0000000
这些结果被认为足够好。外部评级机构提供的信用评级基于定量和定性信息，其中后者最为相关。在我们的案例中，我们仅使用定量公共信息就能预测 68%的评级。
最后，使用rpart.plot包绘制决策树：
prp(DT)
YPCA决策树将如下所示：
模型不是很容易解释吗？
在开始下一部分之前，让我们保存决策树：
save.image("Backup4.RData")
在下一节中，我们将使用另一种有趣的方法，有序逻辑回归，这将能够改进决策树中获得的成果。
有序逻辑回归
正如我们所见，决策树在多分类问题中表现良好。我们还可以遵循其他方法。其中之一是逻辑回归，它对一个问题给出六个可能的结果。然而，这种方法有一些局限性。例如，我们假设目标变量中的类别没有顺序。这意味着目标变量中的不同类别或类是名义的。在评级的情况下，这个假设不一定成立，因为评级分配了一个排名。此外，信用评级之间的差异并不相同，这意味着 AAA 和 AA+评级之间的差异不一定等于 BBB 和 BBB-评级之间的差异。
因此，在这本书的这一部分，我们将实现有序逻辑回归，它假设目标变量中有顺序，并且评级之间的差异不是常数。
该模型可以使用MASS包中的polr函数部署。此函数只需要模型的公式、数据集，在我们的案例中，还需要Hess=TRUE选项。此选项将允许我们计算并可视化模型中变量的标准误差：
library(MASS)
ordered_logistic <- polr(RatingMayT1 ~ ., data = train_macro_trn[,c(variables)], Hess=TRUE)
然后打印出模型的摘要：
summary(ordered_logistic)
## Call:
## polr(formula = RatingMayT1 ~ ., data = train_macro_trn[, c(variables)],
## Hess = TRUE)
##
## Coefficients:
## Value Std. Error t value
## CARA -0.3624 0.2520 -1.4381
## DCPI 0.1432 0.1807 0.7924
## DGDP -0.2225 0.2129 -1.0452
## ILMA 1.5587 0.2592 6.0126
## PSBR 0.6929 0.2209 3.1371
## PUDP -2.8039 0.3886 -7.2145
## TDRA 0.3070 0.2464 1.2461
## YPCA 2.6988 0.7100 3.8011
## MEANWGI 2.2565 0.4707 4.7937
## YearsFromLastDefault 0.8091 0.2191 3.6919
##
## Intercepts:
## Value Std. Error t value
## 1|2 -10.0770 1.1157 -9.0321
## 2|3 -5.6306 0.6134 -9.1789
## 3|4 -2.4390 0.4011 -6.0810
## 4|5 0.4135 0.3615 1.1439
## 5|6 4.8940 0.5963 8.2070
##
## Residual Deviance: 236.9271
## AIC: 266.9271
前面的表格为我们提供了回归系数表。此外，它还显示了不同截距的估计值，有时被称为截点。截距表明预测结果应该在何处被切割，以使数据中观察到的不同信用评级。
此外，该模型为我们提供了残差偏差和AIC指标，这些是用于比较不同模型的有用指标。
在前面的结果中，我们看不到任何表示变量是否显著的p_values，因为在任何回归中通常都不会显示。因此，我们需要计算它们。
p_values可以通过比较t 值与标准正态分布来近似计算。首先，我们将使用以下代码存储我们的系数：
coefs <- coef(summary(ordered_logistic))
print(coefs)
## Value Std. Error t value
## CARA -0.3623788 0.2519888 -1.4380749
## DCPI 0.1432174 0.1807448 0.7923737
## DGDP -0.2225049 0.2128768 -1.0452282
## ILMA 1.5586713 0.2592360 6.0125575
## PSBR 0.6928689 0.2208629 3.1371002
## PUDP -2.8038553 0.3886409 -7.2145133
## TDRA 0.3069968 0.2463570 1.2461463
## YPCA 2.6988066 0.7100112 3.8010760
## MEANWGI 2.2564849 0.4707199 4.7936888
## YearsFromLastDefault 0.8090669 0.2191455 3.6919175
## 1|2 -10.0770197 1.1156894 -9.0321014
## 2|3 -5.6306456 0.6134365 -9.1788566
## 3|4 -2.4389936 0.4010815 -6.0810418
## 4|5 0.4134912 0.3614860 1.1438653
## 5|6 4.8940176 0.5963226 8.2069960
如果我们观察系数的符号，有一些负值。显然，一些变量显示出非直观或意外的符号，但在本例中我们无需对此感到担忧。
模型的系数可能有些难以解释，因为它们是以对数形式缩放的。因此，将之前的原始系数转换为概率比是常见的做法。
概率比将按以下方式获得：
exp(coef(ordered_logistic))
## CARA DCPI DGDP
## 0.69601870 1.15398065 0.80051110
## ILMA PSBR PUDP
## 4.75250240 1.99944358 0.06057607
## TDRA YPCA MEANWGI
## 1.35933662 14.86198455 9.54946258
## YearsFromLastDefault
## 2.24581149
最后，不同变量的p_values被计算并合并到我们获得的系数中：
p_values <- pnorm(abs(coefs[, "t value"]), lower.tail = FALSE) * 2
coefs <- cbind(coefs, "p value" = p_values)
print(coefs)
## Value Std. Error t value p value
## CARA -0.3623788 0.2519888 -1.4380749 1.504128e-01
## DCPI 0.1432174 0.1807448 0.7923737 4.281428e-01
## DGDP -0.2225049 0.2128768 -1.0452282 2.959175e-01
## ILMA 1.5586713 0.2592360 6.0125575 1.826190e-09
## PSBR 0.6928689 0.2208629 3.1371002 1.706278e-03
## PUDP -2.8038553 0.3886409 -7.2145133 5.412723e-13
## TDRA 0.3069968 0.2463570 1.2461463 2.127107e-01
## YPCA 2.6988066 0.7100112 3.8010760 1.440691e-04
## MEANWGI 2.2564849 0.4707199 4.7936888 1.637422e-06
## YearsFromLastDefault 0.8090669 0.2191455 3.6919175 2.225697e-04
## 1|2 -10.0770197 1.1156894 -9.0321014 1.684062e-19
## 2|3 -5.6306456 0.6134365 -9.1788566 4.356928e-20
## 3|4 -2.4389936 0.4010815 -6.0810418 1.194042e-09
## 4|5 0.4134912 0.3614860 1.1438653 2.526795e-01
## 5|6 4.8940176 0.5963226 8.2069960 2.267912e-16
一旦我们开发出我们的模型，我们将预测模型的结果：
Ord_log_pr_train <- cbind(train_macro_trn[,c("CountryName","Year","RatingMayT1")], predict(ordered_logistic, train_macro_trn, type = "probs"))
colnames(Ord_log_pr_train)<-c("Country","year","Observed","X1","X2","X3","X4","X5","X6")
head(Ord_log_pr_train,1)
##Country year Observed X1 X2 X3 X4
1 Austria 2010 6 5.468638e-06 4.608843e-04 0.010757249 0.15316033
## X5 X6
## 1 0.7811701 0.05444599
该模型为每个评级给出了不同的概率。预测评级是预测概率最高的评级。例如，对于 2010 年的奥地利，模型将最高的概率分配给了 5（X5）评级，因此预测评级是 5。
以下代码将预测评级分配给最高概率：
for (j in 1:nrow(Ord_log_pr_train))
{
Ord_log_pr_train$maximaPD[j]<-max(Ord_log_pr_train$X1[j],Ord_log_pr_train$X2[j],Ord_log_pr_train$X3[j],Ord_log_pr_train$X4[j],Ord_log_pr_train$X5[j],Ord_log_pr_train$X6[j])
}
Ord_log_pr_train$Predicted<-ifelse(Ord_log_pr_train$X1==Ord_log_pr_train$maximaPD,1,ifelse(Ord_log_pr_train$X2==Ord_log_pr_train$maximaPD,2,ifelse(Ord_log_pr_train$X3==Ord_log_pr_train$maximaPD,3,ifelse(Ord_log_pr_train$X4==Ord_log_pr_train$maximaPD,4,ifelse(Ord_log_pr_train$X5==Ord_log_pr_train$maximaPD,5,6)))))
让我们看看模型在训练样本中的准确率：
model_assessment(Ord_log_pr_train,"Ordered_logistic")
## notche perc_Ordered_logistic cumulative
## 1 0 0.69047619 0.6904762
## 2 1 0.29761905 0.9880952
## 3 2 0.01190476 1.0000000
如我们所见，该模型能够使用训练样本正确预测 69.05%的信用评级。当使用决策树时，这些结果更好。
让我们看看模型在测试样本中的表现。以下代码给出了每个国家在每个评级水平上的预测概率。预测评级是由概率最高的类别给出的：
Ord_log_pr_test <- cbind(test_macro_trn[,c("CountryName","Year","RatingMayT1")], predict(ordered_logistic, test_macro_trn, type = "probs"))
colnames(Ord_log_pr_test)<-c("Country","year","Observed","X1","X2","X3","X4","X5","X6")
以下代码找到概率最高的评级类别，并将其分配为预测评级：
for (j in 1:nrow(Ord_log_pr_test))
{
Ord_log_pr_test$maximaPD[j]<-max(Ord_log_pr_test$X1[j],Ord_log_pr_test$X2[j],Ord_log_pr_test$X3[j],Ord_log_pr_test$X4[j],Ord_log_pr_test$X5[j],Ord_log_pr_test$X6[j])
}
Ord_log_pr_test$Predicted<-ifelse(Ord_log_pr_test$X1==Ord_log_pr_test$maximaPD,1,ifelse(Ord_log_pr_test$X2==Ord_log_pr_test$maximaPD,2,ifelse(Ord_log_pr_test$X3==Ord_log_pr_test$maximaPD,3,ifelse(Ord_log_pr_test$X4==Ord_log_pr_test$maximaPD,4,ifelse(Ord_log_pr_test$X5==Ord_log_pr_test$maximaPD,5,6)))))
在测试样本中，模型的准确率如下：
model_assessment(Ord_log_pr_test,"Ordered_logistic")
## notche perc_Ordered_logistic cumulative
## 1 0 0.57142857 0.5714286
## 2 1 0.39285714 0.9642857
## 3 2 0.01785714 0.9821429
## 4 3 0.01785714 1.0000000
结果也略逊于决策树模型。
在开始下一节之前，你现在可以保存工作空间：
save.image("Backup5.RData")
在下一节中，我们将使用宏观经济数据来预测国家评级。我们使用的所有变量都是定量变量。在下一节中，我们将使用国家报告来达到相同的目的。
使用欧洲国家报告预测主权评级
根据在“使用宏观经济信息预测国家评级”部分描述的基于宏观经济信息的模型，决策树可以被认为是一种预测主权评级的良好替代方法。
然而，定性信息代表了评级分配中一个重要且低透明度的一部分。在本节中，我们提出了一种仅基于所谓的国家报告的模型，这些报告由欧洲委员会发布。
这些报告主要在二月底发布，包含对欧盟成员国经济和社会挑战的年度分析。
例如，在以下链接中，我们可以下载 2018 年发布的国家报告，ec.europa.eu/info/publications/2018-european-semester-country-reports_en。对于所有 28 个欧盟国家，我们已经从 2011 年到 2018 年下载了它们的国别报告，并将它们转换为文本格式。我们按年份将它们存储在不同的文件夹中，每个文件夹对应一年：
directories <- list.files(path = "../MachineLearning/CountryReports/", pattern = "201", full.names = TRUE)
print(directories)
## [1] "../MachineLearning/CountryReports/2011"
## [2] "../MachineLearning/CountryReports/2012"
## [3] "../MachineLearning/CountryReports/2013"
## [4] "../MachineLearning/CountryReports/2014"
## [5] "../MachineLearning/CountryReports/2015"
## [6] "../MachineLearning/CountryReports/2016"
## [7] "../MachineLearning/CountryReports/2017"
## [8] "../MachineLearning/CountryReports/2018"
让我们创建一个包含每个文件夹中不同报告名称的列表：
txt_files2011<-list.files(path = directories[1], pattern = ".txt", recursive=TRUE,full.names = TRUE)
txt_files2012<-list.files(path = directories[2], pattern = ".txt", recursive=TRUE,full.names = TRUE)
txt_files2013<-list.files(path = directories[3], pattern = ".txt", recursive=TRUE,full.names = TRUE)
txt_files2014<-list.files(path = directories[4], pattern = ".txt", recursive=TRUE,full.names = TRUE)
txt_files2015<-list.files(path = directories[5], pattern = ".txt", recursive=TRUE,full.names = TRUE)
txt_files2016<-list.files(path = directories[6], pattern = ".txt", recursive=TRUE,full.names = TRUE)
txt_files2017<-list.files(path = directories[7], pattern = ".txt", recursive=TRUE,full.names = TRUE)
txt_files2018<-list.files(path = directories[8], pattern = ".txt", recursive=TRUE,full.names = TRUE)
这里，文本文件名存储在一个列表中：
country_reports_list<-do.call(c,list(txt_files2011,txt_files2012,txt_files2013,txt_files2014,txt_files2015,txt_files2016,txt_files2017,txt_files2018))
head(country_reports_list)
## [1] "../MachineLearning/CountryReports/2011/swp_austria_en_0.txt"
## [2] "../MachineLearning/CountryReports/2011/swp_belgium_en_0.txt"
## [3] "../MachineLearning/CountryReports/2011/swp_bulgaria_en_0.txt"
## [4] "../MachineLearning/CountryReports/2011/swp_cyprus_en_0.txt"
## [5] "../MachineLearning/CountryReports/2011/swp_czechrepublic_en_0.txt"
## [6] "../MachineLearning/CountryReports/2011/swp_denmark_en_0.txt"
文件名包含根目录和文件名。让我们尝试将国家名称和报告年份分开：
list<-data.frame(country_reports_list)
list<-data.frame(t(data.frame(strsplit(as.character(list$country_reports_list), "/"))))
list<-list[,(ncol(list)-1):ncol(list)]
row.names(list)<-NULL
list<-cbind(list,country_reports_list)
colnames(list)<-c("Year","file","root")
head(list)
## Year file
## 1 2011 swp_austria_en_0.txt
## 2 2011 swp_belgium_en_0.txt
## 3 2011 swp_bulgaria_en_0.txt
## 4 2011 swp_cyprus_en_0.txt
## 5 2011 swp_czechrepublic_en_0.txt
## 6 2011 swp_denmark_en_0.txt
## root
## 1 ../MachineLearning/CountryReports/2011/swp_austria_en_0.txt
## 2 ../MachineLearning/CountryReports/2011/swp_belgium_en_0.txt
## 3 ../MachineLearning/CountryReports/2011/swp_bulgaria_en_0.txt
## 4 ../MachineLearning/CountryReports/2011/swp_cyprus_en_0.txt
## 5 ../MachineLearning/CountryReports/2011/swp_czechrepublic_en_0.txt
## 6 ../MachineLearning/CountryReports/2011/swp_denmark_en_0.txt
让我们尝试创建一个包含国家名称的列，考虑到每个文件名。例如，如果文件名中包含单词czech，将创建一个新列，其中包含捷克共和国：
list$CountryMapping<-NA
list[grep("austria",list$file),"CountryMapping"]<-"Austria"
list[grep("belgium",list$file),"CountryMapping"]<-"Belgium"
list[grep("bulgaria",list$file),"CountryMapping"]<-"Bulgaria"
list[grep("croatia",list$file),"CountryMapping"]<-"Croatia"
list[grep("cyprus",list$file),"CountryMapping"]<-"Cyprus"
list[grep("czech",list$file),"CountryMapping"]<-"Czech Republic"
list[grep("denmark",list$file),"CountryMapping"]<-"Denmark"
list[grep("estonia",list$file),"CountryMapping"]<-"Estonia"
list[grep("finland",list$file),"CountryMapping"]<-"Finland"
list[grep("france",list$file),"CountryMapping"]<-"France"
list[grep("germany",list$file),"CountryMapping"]<-"Germany"
list[grep("greece",list$file),"CountryMapping"]<-"Greece"
list[grep("hungary",list$file),"CountryMapping"]<-"Hungary"
list[grep("ireland",list$file),"CountryMapping"]<-"Ireland"
list[grep("italy",list$file),"CountryMapping"]<-"Italy"
list[grep("latvia",list$file),"CountryMapping"]<-"Latvia"
list[grep("lithuania",list$file),"CountryMapping"]<-"Lithuania"
list[grep("luxembourg",list$file),"CountryMapping"]<-"Luxembourg"
list[grep("malta",list$file),"CountryMapping"]<-"Malta"
list[grep("netherlands",list$file),"CountryMapping"]<-"Netherlands"
list[grep("poland",list$file),"CountryMapping"]<-"Poland"
list[grep("portugal",list$file),"CountryMapping"]<-"Portugal"
list[grep("romania",list$file),"CountryMapping"]<-"Romania"
list[grep("slovakia",list$file),"CountryMapping"]<-"Slovakia"
list[grep("slovenia",list$file),"CountryMapping"]<-"Slovenia"
list[grep("spain",list$file),"CountryMapping"]<-"Spain"
list[grep("sweden",list$file),"CountryMapping"]<-"Sweden"
list[grep("uk",list$file),"CountryMapping"]<-"United Kingdom"
list[grep("kingdom",list$file),"CountryMapping"]<-"United Kingdom"
list[grep("netherland",list$file),"CountryMapping"]<-"Netherlands"
让我们看看欧盟每个国家有多少个报告：
table(list$CountryMapping)
##
## Austria Belgium Bulgaria Croatia Cyprus
## 8 8 8 6 8
## Czech Republic Denmark Estonia Finland France
## 8 8 8 8 8
## Germany Greece Hungary Ireland Italy
## 8 4 8 8 8
## Latvia Lithuania Luxembourg Malta Netherlands
## 8 8 8 8 8
## Poland Portugal Romania Slovakia Slovenia
## 8 8 8 8 8
## Spain Sweden United Kingdom
## 8 8 8
我们为欧盟的所有国家提供了八个不同的报告，除了克罗地亚（只有6个报告）和希腊（只有4个）。克罗地亚作为欧盟的正式成员国加入是在 2013 年 7 月 1 日。因此，没有 2011 年和 2012 年的报告。至于希腊，2014 年之后没有针对希腊的具体报告。
由于我们打算使用欧洲报告来训练预测信用评级的模型，我们需要选择一些报告来训练模型，其他报告来测试它。在第六章中使用的模型，即在欧洲联盟中可视化经济问题（使用宏观经济信息预测国家评级）部分中使用的相同国家将再次使用。首先，我们需要选择我们之前用于训练模型的那些国家。然后，我们将选定的国家与对应报告所在的位置名称和根目录合并：
train_list<-train_macro[,c("CountryName","Year")]
train_list$year_report<-train_list$Year+1
train_list<-merge(train_list,list,by.x=c("CountryName","year_report"),by.y=c("CountryMapping","Year"),all.x=TRUE)
train_list<-train_list[complete.cases(train_list),]
files_train<-as.vector(train_list$root)
这里是我们将用于训练我们模型的报告示例：
print(head(files_train))
## [1] "../MachineLearning/CountryReports/2011/swp_austria_en_0.txt"
## [2] "../MachineLearning/CountryReports/2012/swd2012_austria_en.txt"
## [3] "../MachineLearning/CountryReports/2013/swd2013_austria_en.txt"
## [4] "../MachineLearning/CountryReports/2014/swd2014_austria_en_0.txt"
## [5] "../MachineLearning/CountryReports/2016/cr2016_austria_en.txt"
## [6] "../MachineLearning/CountryReports/2017/2017-european-semester-country-report-austria-en_1.txt"
同样的程序用于获取验证或测试样本：
test_list<-test_macro[,c("CountryName","Year")]
test_list$year_report<-test_list$Year+1
test_list<-merge(test_list,list,by.x=c("CountryName","year_report"),by.y=c("CountryMapping","Year"),all.x=TRUE)
test_list<-test_list[complete.cases(test_list),]
files_test<-as.vector(test_list$root)
现在我们来看看输出结果：
print(head(files_test))
## [1] "../MachineLearning/CountryReports/2015/cr2015_austria_en.txt"
## [2] "../MachineLearning/CountryReports/2018/2018-european-semester-country- report-austria-en.txt"
## [3] "../MachineLearning/CountryReports/2013/swd2013_belgium_en.txt"
## [4] "../MachineLearning/CountryReports/2011/swp_bulgaria_en_0.txt"
## [5] "../MachineLearning/CountryReports/2013/swd2013_bulgaria_en.txt"
## [6] "../MachineLearning/CountryReports/2014/swd2014_croatia_en.txt"
由于一些国家没有报告，与先前模型中使用的样本大小存在一些差异。以下代码显示了这一点：
print(paste("The number of countries used to train previous model was formed by",nrow(train_macro_trn), "countries",sep=" "))
## [1] "The number of countries used to train previous model was formed by 168 countries"
这是我们将用于训练新模型的国家的数量：
print(paste("The number of countries which we will use to train this new model will be formed by",nrow(train_list), "countries",sep=" "))
## [1] "The number of countries which we will use to train this new model will be formed by 165 countries"
这是用于验证先前模型的国家的数量：
print(paste("The number of countries used to validate previous model was formed by",nrow(test_macro_trn), "countries",sep=" "))
## [1] "The number of countries used to validate previous model was formed by 56 countries"
这是用于训练新模型的国家的数量：
print(paste("The number of countries which we will use to train this new model will be formed by",nrow(test_list), "countries",sep=" "))
## [1] "The number of countries which we will use to train this new model will be formed by 53 countries"
如您所见，差异并不显著。在将文件读入 R 之前，我们将创建一个读取文件的函数。
一旦运行以下函数，我们就可以迭代地读取不同的报告：
Import_txt <- function(txt)
{
x<-as.data.frame(read.delim(txt, header=FALSE, comment.char="#", stringsAsFactors=FALSE))
return(x)
}
将创建两个列表。在每个列表元素中，我们可以找到每个国家的报告：
Reports_train <- lapply(files_train,
function(x)
read.delim(x,
header = FALSE, comment.char="#",
stringsAsFactors = FALSE))
Reports_test <- lapply(files_test,
function(x)
read.delim(x,
header = FALSE, comment.char="#",
stringsAsFactors = FALSE))
在继续之前，可以删除一些不必要的对象，并保存工作区：
rm(list=setdiff(ls(), c("macroeconomic_data","Reports_train","Reports_test","train_list","test_list")))
save.image("Backup6.RData")
报告需要预处理。在提取有用信息或特征以构建我们的模型之前，需要进行数据预处理。
数据清理，或数据预处理，涉及将数据转换为纯文本，然后删除格式、空白、数字、大写字母和停用词。
停用词是指在一种语言中如此常见，以至于它们的信息价值实际上为零的词。由于所有国家报告都可用英文，这些停用词的例子包括介词、限定词和连词。
为了进行这些预处理步骤，加载了tm包，并将报告转换为语料库格式：
library(tm)
docs_train <- as.VCorpus(Reports_train)
docs_test <- as.VCorpus(Reports_test)
创建以下函数以逐个清理我们的报告：
corpus_treatment<-function(corpus)
{
toSpace <- content_transformer(function(x, pattern) {return (gsub(pattern, " ", x))})
corpus <- tm_map(corpus,PlainTextDocument)
corpus <- tm_map(corpus, toSpace, "-")
corpus <- tm_map(corpus, toSpace, ":")
corpus <- tm_map(corpus, removePunctuation)
corpus <- tm_map(corpus, toSpace, "'")
corpus <- tm_map(corpus, toSpace, "'")
corpus <- tm_map(corpus, toSpace, " -")
corpus <- tm_map(corpus,content_transformer(tolower))
corpus <- tm_map(corpus, removeNumbers)
corpus <- tm_map(corpus, removeWords, stopwords("english"))
corpus <- tm_map(corpus, stripWhitespace)
return(corpus)
}
这些报告通过应用以下函数进行转换：
docs_train<-corpus_treatment(docs_train)
docs_test<-corpus_treatment(docs_test)
将进行一项额外的分析，称为词干提取。词干提取过程是指删除词尾以检索词的根（或词干），这在不显著损失信息的情况下减少了数据的复杂性。
因此，动词argue将被缩减为词干argu，无论文本中该词的形式或复杂性如何。因此，其他形式如argued、argues、arguing和argus也将缩减到相同的词干。词干提取过程减少了需要考虑的单词数量，并提供了更好的频率表示。
词干提取过程使用的是SnowballC包：
library(SnowballC)
docs_train <- tm_map(docs_train,stemDocument)
docs_test <- tm_map(docs_test,stemDocument)
在预处理过程之后，考虑国家报告构建了一个矩阵（文档-词矩阵）。这个矩阵的每一行代表每个国家报告，每一列代表在它们上观察到的所有单词。
如果一个词出现在某个国家的报告中，则对应行和列的矩阵条目为 1，否则为 0。当记录文档内的多次出现时，即如果一个词在报告中出现两次，则在相关矩阵条目中记录为 2。
然而，一些提取的单个单词可能缺少原始文本中包含的重要信息，例如词与词之间的依赖关系和高频词周围的上下文。
例如，在报告中提取单词unemployment无法提供足够的信息来解释该术语是积极的还是消极的。因此，从报告中提取的是两个单词的组合，而不是单个单词。
以这种方式，可以在信用度较低的国家中找到一些可能更频繁出现的组合，例如高失业率。
我们将使用名为Rweka的包来提取单词组合：
library(RWeka)
创建以下函数以获取报告中1和2个单词的组合：
options(mc.cores=4)
BigramTokenizer <- function(x) NGramTokenizer(x, Weka_control(min = 1, max = 2))
现在获得了报告的训练和测试样本的文档-词矩阵。只有长度超过3个字符且小于20个字符的单词将被考虑：
tdm2_train <- TermDocumentMatrix(docs_train, control = list(tokenize = BigramTokenizer,wordLengths = c(3,20)))
在验证样本的情况下，我们将考虑在训练样本中找到的单词字典来计算矩阵。测试样本中找到的新单词将不予考虑。请记住，我们不能使用测试样本中找到的单词来开发模型，这意味着，现在我们应该假装测试样本不存在来训练算法。
因此，当在报告的测试样本上计算文档-术语矩阵时，我们应该在我们的函数中添加一个新参数：dictionary = Terms(tdm2_train)。这在上面的代码中显示：
tdm2_test <- TermDocumentMatrix(docs_test, control = list(dictionary = Terms(tdm2_train),tokenize = BigramTokenizer,wordLengths = c(3,20)))
让我们分析结果矩阵：
tdm2_train
## <<TermDocumentMatrix (terms: 609929, documents: 165)>>
## Non-/sparse entries: 2034690/98603595
## Sparsity : 98%
## Maximal term length: 20
## Weighting : term frequency (tf)
tdm2_test
## <<TermDocumentMatrix (terms: 609929, documents: 53)>>
## Non-/sparse entries: 504543/31821694
## Sparsity : 98%
## Maximal term length: 20
## Weighting : term frequency (tf)
第一行表示每个样本中不同术语的数量和报告的数量。每个矩阵的列数或术语数是相同的，属于我们在训练样本中找到的总单词列表。
总共有 609,929 个术语至少在国家报告的训练样本中出现过一次。
此外，在训练矩阵中，有 98,603,595 个单元格的频率为 0，而 2,034,690 个单元格具有非零值。因此，所有单元格中有 98%是零。这种情况在文本挖掘问题中非常常见。当一个矩阵包含许多零值单元格时，它被认为是稀疏矩阵。
removeSparseTerms函数将删除使用频率较低的单词，只留下语料库中最常用的单词。
在我们的情况下，我们将矩阵减少到最多保留 75%的空空间。这个过程必须只应用于我们稍后用于训练模型的那些数据：
tdm2_train2 <- removeSparseTerms(tdm2_train, 0.75)
tdm2_train2
## <<TermDocumentMatrix (terms: 6480, documents: 165)>>
## Non-/sparse entries: 589204/479996
## Sparsity : 45%
## Maximal term length: 20
## Weighting : term frequency (tf)
我们现在用于训练模型的矩阵有 6,480 个术语。
让我们观察我们的矩阵看起来是什么样子：
print(as.matrix(tdm2_train2[1:10, 1:4]))
## Docs
## Terms character(0) character(0) character(0) character(0)
## gj 0 0 0 0
## mwh 0 0 0 0
## ± ± 0 0 0 0
## à vis 1 1 1 5
## abil 0 1 1 2
## abl 0 1 1 0
## abl afford 0 0 0 0
## abolish 1 1 3 1
## abroad 2 4 3 0
## absenc 0 0 0 0
例如，单词abroad在第一份报告中出现了2次，在第二份报告中出现了4次。只需记住，在之前的步骤中已经进行了词干提取，所以只显示单词的根。矩阵中也包含单个单词和两个单词的组合。
在前面代码中显示的报告名称是按照我们最初用来导入报告的列表顺序排列的。具体来说，前四个文档属于以下：
print(head(as.character(train_list$root),4))
## [1] "../MachineLearning/CountryReports/2011/swp_austria_en_0.txt"
## [2] "../MachineLearning/CountryReports/2012/swd2012_austria_en.txt"
## [3] "../MachineLearning/CountryReports/2013/swd2013_austria_en.txt"
## [4] "../MachineLearning/CountryReports/2014/swd2014_austria_en_0.txt"
也可能获得术语及其频率的完整列表：
freq <- rowSums(as.matrix(tdm2_train2))
ord <- order(freq,decreasing=TRUE)
这是出现频率最高的术语列表：
freq[head(ord,20)]
## gdp rate increas market labour sector tax growth public
## 20566 16965 15795 15759 15751 15381 14582 14545 14515
## year employ energi invest measur govern debt bank term
## 13118 12481 12330 12027 11341 10903 10854 10668 10470
## averag social
## 10059 10051
为了可视化，我们可以创建一个包含我们文档中最频繁单词的词云图。为此，我们可以使用wordcloud包：
library(wordcloud)
set.seed(1234)
wordcloud(row.names(tdm2_train2), freq = freq, max.words=200,min.freq=4000,scale=c(2,.4),
random.order = FALSE,rot.per=.5,vfont=c("sans serif","plain"),colors=palette())
结果看起来就是这样：
看起来很漂亮，对吧？
现在，我们需要将之前的矩阵转换为以国家为行、术语为列。简而言之，我们需要转置训练和测试矩阵：
tdm2_train2 <- as.matrix(tdm2_train2)
dim(tdm2_train2)
## [1] 6480 165
tdm2_train2 <- t(tdm2_train2)
tdm2_test2<-as.matrix(tdm2_test)
tdm2_test2 <- t(tdm2_test2)
rm(tdm2_test)
再次，一些无关的对象被移除了：
rm(list=setdiff(ls(), c("macroeconomic_data","train_list","test_list","tdm2_test2","tdm2_train2")))
并且工作区再次保存为备份：
save.image("Backup7.RData")
到目前为止，我们已经处理了我们的报告，并从中提取了一些特征、术语和单词组合。尽管如此，目标变量，即国家评级，并不在我们新的数据集中。信用评级仅在macroeconomic_data样本中存在。
在下面的代码中，我们将为最近创建的训练和验证矩阵添加信用评级：
train_list<-merge(train_list[,c("Year","CountryName","year_report")],macroeconomic_data[,c("CountryName","Year","RatingMayT1")],by=c("CountryName","Year"),all.x=TRUE)
test_list<-merge(test_list[,c("Year","CountryName","year_report")],macroeconomic_data[,c("CountryName","Year","RatingMayT1")],by=c("CountryName","Year"),all.x=TRUE)
training <- cbind(train_list,tdm2_train2)
validation <- cbind(test_list,tdm2_test2)
由于我们模型中要训练的特征数量相当高（超过 6,000 个），我们将评估我们的特征与信用评级的相关性，以帮助排除其中的一些。
首先，我们将创建一个包含我们的术语列表和与信用评级相关性的数据框。前三个变量必须排除。以下代码显示了这一点：
head(colnames(training),7)
## [1] "CountryName" "Year" "year_report" "RatingMayT1" " gj"
## [6] " mwh" "± ±"
现在我们有了相关性，让我们按降序排列它们：
correlations<-data.frame(correlations)
colnames(correlations)<-c("word","correlation")
correlations$abs_corr<-abs(as.numeric(as.character(correlations$correlation)))
correlations<-correlations[order(correlations$abs_corr,decreasing = TRUE),]
correlations = matrix("NA",nrow=(ncol(training)-4),2)
ncolumns<-ncol(training)
for (i in 5:ncolumns)
{
correlations[i-4,1]<-colnames(training[i])
correlations[i-4,2]<- as.numeric(cor(training[,i],as.numeric(as.character(training[,"RatingMayT1"]))))
}
这里是信用评级相关性最高的前 10 个变量：
head(correlations,10)
## word correlation abs_corr
## 3245 judici -0.495216233392176 0.4952162
## 1175 court -0.4939081009835 0.4939081
## 132 administr -0.470760214895828 0.4707602
## 3710 migrant 0.460837714113155 0.4608377
## 1343 delay -0.46038844705712 0.4603884
## 468 background 0.455839970556903 0.4558400
## 116 adequ -0.445062248908142 0.4450622
## 2811 immigr 0.428818668867468 0.4288187
## 3246 judici system -0.42745138771952 0.4274514
## 6106 undeclar -0.419206156830568 0.4192062
显然，来自诸如司法等单词的judici词根与信用评级高度相关。负号表示在报告中出现频率非常高的特定单词的国家信用质量较低。
我们将只使用前 1,000 个单词来训练我们的模型。这里创建了一个包含前 1,000 个术语的列表：
list_vars<-dput(as.vector(correlations$word[1:1000]))
在训练模型之前，让我们再次保存工作空间：
save.image("Backup8.RData")
是时候训练模型了。选定的模型是一个纯 Lasso 模型，因为已经证明这种模型在列数或特征数较多的情况下效果良好，它作为一种变量选择的方法。
这种方法已经在第五章：预测银行失败 - 多变量分析中使用过，使用了h2o包。这次，我们仅为了学术目的使用glmnet包，目的是让读者可以应用不同的解决方案：
library(glmnet)
glmnet 包需要一个包含变量的矩阵和一个包含类别标签或目标值的向量。
让我们确保我们的目标变量是一个factor：
training$RatingMayT1<-as.factor(training$RatingMayT1)
validation$RatingMayT1<-as.factor(validation$RatingMayT1)
依赖变量和独立变量存储在不同的对象中，如下面的代码所示，以训练模型：
xtrain<-training[,list_vars]
ytrain<-training$RatingMayT1
与前面的代码一样，相同的步骤在验证样本中执行：
validation$RatingMayT1<-as.factor(validation$RatingMayT1)
xtest<-validation[,list_vars]
ytest<-validation$RatingMayT1
我们还将使用cv.glmnet函数在训练过程中，该函数自动执行网格搜索以找到 Lasso 算法中所需的λ的最佳值。
此函数中最重要的参数如下：
-
y：我们的目标变量，在本例中，是信用评级。 -
x：一个包含我们特征所有独立变量的矩阵。 -
alpha：在我们的情况下，值为1表示模型是 Lasso。 -
family：我们的响应变量的类型。如果目标变量只有两个水平，则应将family定义为binomial。在我们的情况下，由于我们的目标变量显示超过两个水平，因此应将family指定为multinomial。 -
type.multinomial：如果grouped，则对变量的multinomial系数使用分组 Lasso 惩罚。默认为ungrouped。 -
parallel：如果TRUE，算法将以并行方式处理。这意味着算法将不同的任务分割并同时执行，显著减少训练时间。
下面是使用当前数据的此函数的应用：
set.seed(1234)
ModelLasso <- cv.glmnet(y = ytrain, x=data.matrix(xtrain[,list_vars]), alpha=1,family='multinomial',type.multinomial = "grouped",parallel=TRUE)
在执行此代码的过程中，出现了一个警告信息：one multinomial or binomial class has fewer than 8 observations; dangerous ground。
问题在于我们对于目标变量中的所有类别都没有足够的观测数据。我们可以通过运行以下代码来检查目标变量中不同类别的数量：
table(ytrain)
## ytrain
## 1 2 3 4 5 6
## 5 23 33 32 34 38
对于评级1，只有5个观测值。因此，对于这个类别，可能不会期望有任何稳定的估计。
一种可能的解决方案是将评级1和2合并到同一个评级类别中：
ytrain<-gsub("1","2",ytrain)
ytest<-gsub("1","2",ytest)
现在，问题应该不会出现了：
table(ytrain)
## ytrain
## 2 3 4 5 6
## 28 33 32 34 38
set.seed(1234)
ModelLasso <- cv.glmnet(y = ytrain, x=data.matrix(xtrain[,list_vars]), alpha=1,family='multinomial',type.multinomial = "grouped")
模型训练完成后，以下图表有助于找到减少模型误差的lambda参数：
plot(ModelLasso)
根据以下图表，最优对数值大约为-3：
可以通过检查代码中的lambda_min变量来查看确切值：
log(ModelLasso$lambda.min)
## [1] -3.836699
正则化方法的目标是在准确性和简单性之间找到一个平衡点，这意味着要获得一个具有最小系数数量且也能给出良好准确率的模型。在这方面，cv.glmnet函数也有助于找到误差在最小误差一倍标准差内的模型。
这个lambda值可以在lambda.1se变量中找到。这个值将被选为我们模型的最终lambda值：
best_lambda <- ModelLasso$lambda.1se
print(best_lambda)
## [1] 0.05727767
现在，是时候评估我们模型的准确率了。首先，让我们看看训练样本：
predictLASSO_train <- predict(ModelLasso, newx = data.matrix(xtrain[,list_vars]),
type = "class", s = ModelLasso$lambda.1se)
predictLASSO_train<-as.data.frame(cbind(training[,1:2],ytrain ,predictLASSO_train))
colnames(predictLASSO_train)<-c("Country","Year","Rating","Prediction")
以下表格是训练样本的结果表：
table(predictLASSO_train$Rating,predictLASSO_train$Prediction)
## 2 3 4 5 6
## 2 27 0 0 0 1
## 3 1 32 0 0 0
## 4 0 1 30 0 1
## 5 0 0 0 33 1
## 6 0 0 0 1 37
现在，让我们看看验证样本的准确率：
predictLASSO_test <- predict(ModelLasso, newx = data.matrix(xtest),
type = "class", s = ModelLasso$lambda.1se)
predictLASSO_test<-as.data.frame(cbind(validation[,1:2],ytest ,predictLASSO_test))
colnames(predictLASSO_test)<-c("Country","Year","Rating","Prediction")
以下表格是验证样本的结果表：
table(predictLASSO_test$Rating,predictLASSO_test$Prediction)
## 2 3 4 5 6
## 2 5 3 1 0 1
## 3 1 7 1 0 0
## 4 0 0 7 0 3
## 5 0 1 1 8 2
## 6 0 0 0 2 10
考虑到我们使用了国家报告，结果似乎已经足够好。正如我们使用宏观经济数据训练模型时做的那样，我们将使用以下函数计算正确分类国家的百分比：
model_assessment<-function(data,model)
{
data$Observed<-as.numeric(as.character(data$Rating))
data$Predicted<-as.numeric(as.character(data$Prediction))
data$df<-abs(as.numeric(data$Predicted)-as.numeric(data$Observed))
comparison<-as.data.frame(table(data$df))
comparison$perc<-comparison$Freq/nrow(data)
colnames(comparison)<-c("notch","N",paste("perc_",model,sep=''))
comparison$N<-NULL
return(comparison)
}
让我们运行这个模型的评估：
model_assessment(predictLASSO_train,"Train_LASSO")
## notch perc_Train_LASSO
## 1 0 0.963636364
## 2 1 0.024242424
## 3 2 0.006060606
## 4 4 0.006060606
model_assessment(predictLASSO_test,"Test_LASSO")
## notch perc_Test_LASSO
## 1 0 0.69811321
## 2 1 0.18867925
## 3 2 0.09433962
## 4 4 0.01886792
Lasso 模型在验证样本中能够正确预测 69.81%的国家。由此得到的模型在仅使用宏观经济数据获得的结果上略有改进，达到了 67.86%的准确率。
最后，评估国家报告中出现并决定国家信用评级的最重要的术语是非常有趣的。
以下函数用于提取模型的系数：
coefs<-coef(ModelLasso, s = "lambda.1se")
结果是一个列表，列出了每个评级级别的不同系数。例如，信用评级 1 和 2（这些类别在此部分之前已合并）的系数被获得。这将在以下代码中显示：
coefs2<-coefs$`2`
list_coefs2<-as.data.frame(coefs2@Dimnames)
colnames(list_coefs2)<-c("variable","id")
list_coefs2$id<-as.numeric(row.names(list_coefs2))-1
aux_coefs2<-cbind(as.data.frame(coefs2@i),as.data.frame(coefs2@x))
colnames(aux_coefs2)<-c("id","coefficient")
list_coefs2<-merge(list_coefs2,aux_coefs2,by.x="id")
rm(coefs2,aux_coefs2)
这里显示了一些最相关的术语：
head(list_coefs2[order(list_coefs2$coefficient,decreasing = TRUE),],10)
## id variable coefficient
## 18 69 financ need 0.24991828
## 37 192 personnel 0.13635379
## 44 305 outflow 0.11243899
## 15 51 energi sector 0.06854058
## 24 97 minimum incom 0.05821313
## 39 216 gross extern 0.05237113
## 10 37 resolut 0.04807981
## 72 700 analyt 0.03036531
## 75 774 healthcar sector 0.02997181
## 26 102 social benefit 0.02572995
正面的迹象表明，一个术语在国别报告中出现的频率越高，该国的信用质量就越低。
让我们检查这是否被正确观察到。例如，模型检查了 2018 年塞浦路斯国别报告中包含financ need的一些句子。以下是报告的三个部分：
-
塞浦路斯似乎不面临立即的财政压力风险，这主要得益于其有利的财政地位。这主要归功于一般政府财政平衡和初级平衡的改善、低融资需求以及相对较低短期一般政府债务。这些因素超过了仍然相当大的公共债务。然而，宏观金融方面的短期风险仍然显著。
-
融资渠道有限和降低债务的需求仍然抑制了私营部门的投资。
-
公共债务显著下降，但仍然很高，2017 年约为 GDP 的 99%。高公共债务使塞浦路斯容易受到金融或经济冲击。然而，在经济调整计划期间，外部债权人提供的长期低息债务的大比例、当前低主权债券收益率以及相对较低的中期融资需求减轻了再融资风险。
在这三个部分中，去除了停用词，这也是找到financ need的原因。
对于最佳的评级类别，也可以得到不同的系数。这可以通过以下代码实现：
coefs6<-coefs$`6`
list_coefs6<-as.data.frame(coefs6@Dimnames)
colnames(list_coefs6)<-c("variable","id")
list_coefs6$id<-as.numeric(row.names(list_coefs6))-1
aux_coefs6<-cbind(as.data.frame(coefs6@i),as.data.frame(coefs6@x))
colnames(aux_coefs6)<-c("id","coefficient")
list_coefs6<-merge(list_coefs6,aux_coefs6,by.x="id")
rm(coefs6,aux_coefs6)
这里是我们找到的一些最相关的术语：
head(list_coefs6[order(list_coefs6$coefficient,decreasing = TRUE),],10)
## id variable coefficient
## 45 309 remaind 0.22800169
## 1 0 (Intercept) 0.20122381
## 7 20 govern balanc 0.15410796
## 81 899 stimulus 0.11734883
## 82 918 europ strategi 0.06968609
## 17 57 interest payment 0.05516403
## 49 367 fiscal posit 0.04272709
## 65 568 contribut rate 0.03101503
## 38 207 decad 0.03063200
## 2 6 background 0.03029957
还获得了一些示例报告的句子。例如，对于 2018 年的德国，以下句子包含govern balanc的组合：
-
德国一直改善其政府平衡，从 2014 年开始转变为盈余。
-
积极的政府平衡也反映在政府债务的下降上，2015 年达到 70.9%，进一步下降到 2016 年的 68.1%。
最后，为了以防万一你想以后使用，请备份所有你的模型：
save.image("Backup9.RData")
摘要
在本章中，你学习了文本挖掘和主题提取的一些基本概念。你现在应该知道如何读取文本文件并处理原始文本以获取有用的常用词。此外，你现在能够在自己的问题中使用以文本格式收集的信息。
根据您想要解决的数据量和问题类型，您现在可以应用本书中先前使用过的各种技术，无论是简单的还是复杂的。
最后，考虑到本章内容，你已经准备好深入探索其他更近和更有前景的技术，例如 word2vec 和 doc2vec，这两种都是高级技术，允许你在文本和文档中发现相关信息或主题。如果你对此感兴趣，可以进一步研究这些主题。
我希望你对机器学习有了深入的了解，并且这本书帮助你开始了使用机器学习解决问题的旅程。感谢阅读，祝您一切顺利！
