《R 语言》实验报告

PAGEPAGE21《R语言》实验报告专业班级： XXXX 成员姓名：XXXXXXXXXXXXX 指导教师： XXX 日期： XXXX-XX-XX 一、问题重述 31.1问题背景 3二、问题分析 3问题一的分析 3问题二的分析 3三、数据预处理 4数据选取 4缺失值处理 5数据统一化 5数据去极值 7数据标准化 7四、模型的建立与求解 7问题一模型的建立与求解 7模型的建立 8模型的求解 9模型的结果 10问题二模型的建立与求解 模型的建立 模型的求解 13模型的结果 14五、模型评价与改进 15六、附录 16附录1 16附录2 18一、问题重述问题背景市场信息中提取出有效指标，制订交易策略，是一个具有挑战性的工作。问题提出202171420221285据信息，解决以下问题：问题一：在所提供的49项指标中，筛选出与“数字经济”板块有关的主要指标。问题二：建立模型对每5分钟的“数字经济”板块指数进行预测。二、问题分析问题一的分析49较强的指标后进行主成分分析，得到最终的指标。问题二的分析5开盘价与收盘价之差的和定义为差价，同时将差价以零为分界线划分为两部分，采用LogisticLogistic过程复杂。Lasso0LassoZou2006AdaptiveLassoLassoLogisticLogistic在的指标多重共线性和计算复杂等问题。三、数据预处理数据选取本文选取“数字经济（CSI:931582）”49个股票量化指标作为2021714202212849个量化指标分为五大类：宏观市场指标、国内股票市场指标、技术指标、国际股票市场指标、（见iFinDEXCEL、SPSS与R。1宏观市场指标采购经理指数、社会消费品零售总额、居民消费价格指数、人民币贷款利率国内股票市场指标上证综合指数成交量、上证综合指数成交金额、沪市股票流通市30050050A深证综合指数、股票市场总值技术指标VMA、VMACD、ARBR、OBV、BBI、DMA、MA、EXPMA、MTM、MACD、BIAS、KDJ、RSI、BOLL国际股票市场指标道琼斯工业指数、纳斯达克综合指数、标准普尔500指数、美国225100CAC40AEXMIB指数、美元/人民币汇率、欧元/美元汇率其他板块信息数字媒体、数字孪生、快手概念、互联网电商、互联网缺失值处理记{i}nn表示量,i和mA=(1,…,n,𝑌)A的每一行为一个观即利用补齐的办法，构造一个能够充分反应原始数据，且不存在数据缺失的测量矩阵。经检查（见2），680处数据且缺失指标数据量较多，故采用分段插值的方式对其进行补全。数据统一化5以日为间隔的，故需将数据的标度统一为日。y=0.0257x3-5.4752x2+253.48x+38628R2=0.94335y=0.0257x3-5.4752x2+253.48x+38628R2=0.9433430004100039000370003500033000

0 20 40 60 80 100 120 140 160图1“社会消费品零售总额”数据拟合2指标名称缺失值社会消费品零售总额22/1/31（39130.1）、21/2/28（35626.7）、21/1/31（37650.9）居民消费价格指数22/1/31（100.9）人民币贷款利率22/1/31（7.5）股票市场总值22/1/31（839226）道琼斯工业平均指数、21/10/31（35866.6）、21/9/6（35234.5）、21/7/31（34886.8）纳斯达克综合指数、21/10/31（15547.1）、21/9/6（15368.9）、21/7/31（14676.8）标准普尔500指数22/1/17（4619.9）、21/12/24（4758.4）、21/11/25（4648.0）、21/10/31（4609.5）、21/9/6（4527.7）、21/7/31（4391.2）美国证交所指数22/1/17（3645.7）、21/12/24（3459.5）、21/11/25（3452.4）、21/10/31（3453.0）、21/9/6（3071.5）、21/7/31（3015.5）香港恒生指数、、21/7/31（26098.4）东京日经225指数、、21/9/23（29944.1）、21/9/20（30169.8）、21/8/9（27854.0）、21/7/31（27532.3）、21/7/23（27738.1）、21/7/22（27643.0）伦敦金融时报100指数22/1/3（7444.8）、21/12/28（7404.4）、21/12/27（7388.2）、21/10/31（7263.0）、21/8/30（7133.8）、21/7/31（7057.0）法国巴黎CAC40指数21/10/31（6861.8）、21/7/31（6644.3）荷兰AEX指数21/10/31（811.3）、21/7/31（756.2）俄罗斯RTS指数21/10/31（1855.4）、21/7/31（1630.8）意大利MIB指数、21/7/31（25364.5）美元/人民币汇率、、21/7/31（6.4）欧元/美元汇率、、21/7/31（1.1）数据去极值（3𝜎准则限的极端值用上下限值代替。记

=fi

mi} (𝑖=1,…,𝑛)为所考虑量化指标的数据，各指标数据为𝑎=i=1(i1,…,imim1,…,m̅i=i−̅(𝑖=1,…,𝑚，从而得到标准差为1𝜎=𝑚

m(i−̅)2i=1某个测量值𝑥i的剩余误差𝑣i满足|𝑣i|=|𝑥i−𝑥̅|>3𝜎，则认为𝑥i有比较大的误差，将应予代替。数据标准化01𝑋=(1,…,m为原始数据，𝑋∗=(𝑥∗,…,𝑥∗为标准化的数据，转化公式为1 m𝑋∗=i−𝜇i(𝑖=1,…𝑛)i i𝜇iiii个指标数据的标准差。四、模型的建立与求解问题一模型的建立与求解ICICIC（informationcoefficient）是一个用于衡量预测值优劣的绩效指维因子数据的主因子。模型的建立两个变量之间的皮尔逊相关系数定义为两个变量之间的协方差和标准差的商：𝑐𝑜𝑣(𝑋,𝑌) 𝐸(𝑋𝑌)−𝐸(𝑋)𝐸(𝑌)X,y=

=y

𝜎X

y上式定义了总体相关系数，常用希腊小写字母𝜌作为代表符号。估算样本的协方差和标准差，可得到皮尔逊相关系数，常用英文小写字母r代表：∑m(Xi−)(i−)r= i=1 ✓∑mi=1✓∑

(Xi−)2

(i−)2r亦可由(Xi,Yi)样本点的标准分数均值估计，得到与上式等价的表达式：1 m X−X

Y−Yr=m−

I（iσXi=1

）（i ）σYXi--、及σX样本的标准分数、样本平均值和样本标准差。ax X i引起该系数的改变，即它该变化的不变量(由符号确定)。也就是说，我们如果把𝑋移动到𝑎𝑏𝑋和把𝑌移动到𝑐𝑑𝑌，a、b、cd是常数，并不会改变两个变量的相关（该结论在总体和样本皮尔逊相关系数中都成立会改变相关系数：∑ii−𝑚̅- 𝑚∑ii−∑i∑ixy=

(𝑚−

y

= 𝑚∑𝑥2−(∑𝑥)2𝑚∑𝑦2−(∑𝑦)2i i i inkknk维特征。记x=(x1xm)m维样本数据，样本均值与方差分别为m1= I

𝑆2=

mI(𝑥

−̅)2样本𝑥y

𝑚 ii=1

𝑚−

ii=1m1cv(x,)=m−1I(xi−-)(i−-)i=1协方差为正则正相关，为负则负相关，为零则不相关。设XX1Xn)n维随机变量，Xim维样本数据，则Xmn列的矩阵。Xn阶的方阵cov(X1,X1) cov(X1,X2)cov(X)=cov(X2,X1) cov(X2,X2)⋮ ⋮cov(Xn,X1) cov(Xn,X2)

… cov(X1,Xn)… cov(X2,Xn)⋱ ⋮… cov(Xn,Xn)设cv(X)的特征值按从大到小排列为1≥···≥𝜆n，其对应的单位特征向量分别为1,...,niiii=(X1,…,Xn)i前k个主成分的累计方差贡献率记为E(k)=

k∑∑λλj=1j∑∑λλnj=1j累积方差贡献率越高，相应的主成分能够涵盖的信息量就越大。本文确定k使得E(k)>90%，同时要求相应的特征值大于1。模型的求解SPSS4914（表4300MTMBIASRSI与因变量具有较强的相关性。KMOBartlett（见3），KMO，KMO0.736，BartlettP0.01，设，各变量间具有相关性，主成分分析有效，程度为一般。表3KMO检验和Bartlett的检验KMO值0.736Bartlett球形度检验近似卡方972.163df15.000p0.000***4差价沪深300指数 0.259(0.003***)上证综合指数 0.161(0.068*)创业板指数 0.285(0.001***)上证50指数 0.180(0.041**)上证A股指数 0.161(0.067*)深证成份指数 0.277(0.001***)深证综合指数 0.213(0.015**)ARBR 0.179(0.042**)MTM 0.333(0.000***)BIAS 0.458(0.000***)KDJ 0.163(0.064*)RSI 0.497(0.000***)纳斯达克综合指数 0.181(0.039**)东京日经225指数 0.152(0.085*)模型的结果SPSS主成分分析后，得到如下5（主成分载荷出成分得分，得出因子公式，其计算公式为：线性组合系数*（方差解释率/释率），最后将其归一化即为因子权重得分。通过分析成分矩阵（见表6），1290.556%。模型的公式为：F1=0.184×沪深300指数+0.194×创业板指数+0.197×深证成份指数+0.196×MTM+0.212×BIAS+0.206×RSIF2=0.395×沪深300指数+0.328×创业板指数+0.426×深证成份指数-0.385×MTM-0.343×BIAS-0.351×RSIF=(0.706/0.906)×F1+(0.2/0.906)×F25成分初始特征根提取平方和载入特征根方差百分比累积特征根方差百分比累积14.23570.579%70.579%4.23570.579%70.579%21.19919.976%90.556%1.19919.976%90.556%30.3295.488%96.043%40.1442.404%98.447%50.0641.066%99.513%60.0290.487%100.0%6名称成分成分1成分2沪深300指数0.1840.395创业板指数0.1940.328深证成份指数0.1970.426MTM0.196-0.385BIAS0.212-0.343RSI0.206-0.351问题二模型的建立与求解模型的建立“亏损和“盈利01设(xi,yi)为m个独立的观测样本值，xi表示第i天的n个指标的观测值，yi为响应变量，i∈{0,}i=𝑝(i=1iiLogit变化，得：log𝑖𝑡𝜋

=ln(

m)=𝛽+I𝑥T𝛽=

n+I𝑥𝛽i 1−

。 i=1

。 i=1。i=(i1,i2,…,in)Tn𝛽=(1,2,…,n)TLogistic过最大化对数似然值来估计参数。则上式的对数似然函数为：m m𝑙(𝛽,𝛽)=I{𝑦ln(𝜋)+(1−𝑦)ln(1−𝜋))}=I𝑦(𝛽

+𝑥T𝛽)−log1+𝑒{o+xT{J}。 i ii=1

i i i 。 i ii=1将上式对𝛽求导，令导数等于0，可得：𝜕𝑙

m=−

i

𝑒{o+{ii

m+I𝑦𝑥 =0𝜕𝛽

i=1

1+𝑒oi

ii进而有𝑃(𝑥;𝛽;𝛽)=oi

mI(i−𝑝(i;。;i))i=0i=1。i 。

1+{o+{ixi由于上式无法求解析解，故选用经典的梯度下降算法。按照经典梯度下降算法，定义似然函数为𝑓(𝛽)，更新公式为：𝛽+)=𝛽)−𝐻-1(𝛽())𝑓(𝛽())其中∇𝑓f的梯度，它的偏导数矢量为∇𝑓叮,叮]。ＨfHessian矩阵，{J1

{Jk它的二阶偏导数矩阵为𝐻=2叮{i{j

。容易看出，当k＝２时，即为二分类模型。综上，可以得到第ｉ个样本被分到ｙ＝１的概率估计值为：i＾=i

xp(。+

i)=11+x(。=1

n++

i

𝛽)若＾i>0.。LogisticLogistic（PenalizedLogisticLogisticL1到如下式子：mPLR=−I{iln(i)+(1−i)ln(1−i))}+λP(β)i=1β估计通过最小化上式可得：𝛽^

m=𝑎𝑟𝑔min−{𝑦ln(𝜋)+(1−𝑦)ln(1−𝜋))}+λP(β){J i i i ii=1其中λP(β)为惩罚项，由正的调节参数λ决定，λ控制着模型拟合和惩罚项的影响之间的权衡，从而避免训练数据中过度拟合的问题。我们不惩罚截距项，若惩罚项为λP(β)=λ∑pIβI，则可得系数的Lasso估计为:j=1 j𝛽^

n=𝑎𝑟𝑔min−

pln(𝜋)+(1−𝑦)ln(1−𝜋))}+λIβI{J i i=1

i i jj=1当λ＝0MLE估计，λλ→∞，0差之间的正确平衡，以便最小化误分类错误。所以选取适当的λ是非常重要的，通常使AIC，BIC，广义交叉验证。模型的求解首先，根据差价将因变量划分为“盈利”和“亏损”，对因变量进行频数统计得到表7。由于二者比例接近1：1，故模型受选项影响较小。7因变量差价选项盈利频数69百分比53.08%亏损6146.92%总计130100%AICBIC值。第一个模型以第一问筛LogisticLassoLogistic的结果8效果，但数据却造成了损失。在第三个模型中，Lasso方法在自变量中剔除了主成分1使得模型的解释性更好。8变量名称模型选择模型1模型2模型3x1沪深300指数0.31890.2950x2创业板指数-0.2581-0.2576x3深证成份指数-0.1809-0.1639x4MTM-1.2956-1.2969x5BIAS-0.2154-0.1257x6RSI2.37842.3022x7主成分10.2146--x8主成分2-0.51930.01451截距项0.19220.13110.1906AIC161.2494172.4127BIC181.3221181.0153模型的结果9139828131型的解释性更好。9模型1模型2模型3数量988198平方误差168.3267113.4467165.2959五、模型评价与改进AdaptiveLassoLogistic回归在惩罚项上加权重，有Oracle性质n p𝛽^

=𝑎𝑟𝑔min−I{𝑦ln(𝜋)+(1−𝑦)ln(1−𝜋))}+λIωIβI{J i i=1

i i j jj=1其中ω=(ω1ω2ωp)T是ｐ×1jωj≥0Lasso给大的系数，即大的系数接受小的惩罚，小的系数接受大的惩罚。j=(IjI)yγ≥0^λ和γ两个调节参数，为了计算方便，令γ=1λ即可。AdaptiveLassoλOracleLogistic提出的最小二乘近似方法。n(𝛽)是负的对数似然函数，假设n(𝛽)关于𝛽𝛽是通过最小化n(𝛽矿aylor展开式为：-1

-1

1 T1𝑛 n(𝛽)≈

n(𝛽)+

n(𝛽)

(𝛽−𝛽)+ (𝛽−𝛽2

[𝑛n(𝛽)(𝛽−𝛽)其中̇()̈()分别是n()矿n(𝛽)̇()=0，则上式可以简化成-1

1 T1𝑛 n(𝛽)≈

n(𝛽)+2(𝛽−𝛽)[𝑛n(𝛽)(𝛽−𝛽)忽略常数项

()1，进一步简化为2-1 T1𝑛 n(𝛽)≈(𝛽−𝛽)[𝑛n(𝛽)(𝛽−𝛽)特别地E{1𝐿̈()}=∑-1（∑为的渐近协方差阵），所以可以＾-1 1 作n n =nn(𝛽)为∑-1(𝛽−)T＾-1(𝛽−)1𝐿̈()的一nn个简单近似。最终得到统一的基于最小二乘近似的AdaptiveLassop^=𝑎𝑟𝑔min(𝛽−)T＾1(𝛽−)+𝜆I矿jIβjI{Jj=1六、附录附录1y=-0.0018x2y=-0.0018x2+0.206x+48.815R2=0.56825550450 20 40 60 80 100 120 140 160“采购经理指数”数据拟合y=9E-06x3y=9E-06x3-0.0018x2+0.0902x+100.69R2=0.8188102.5102101.5101100.50 20 40 60 80 100 120 140 160“居民消费价格指数”数据拟合8.58.38.17.9

y=7E-08x4-2E-05x3+0.0023x2-0.0693x+7.5825R²=0.6077.77.57.37.16.96.76.50 20 40 60 80 100 120 140 160“人民币贷款利率”数据拟合LassoLogistic模型的系数变化图2021812日K附录2#绘制K线图 3. library(tibble)不可缺少3. library(tibble)不可缺少4. library(magrittr)5. 5. library(quantmod)7. shares7. sharesread.csv("F:/shares_data.csv"8. shares$timeas.Date(shares$time9. 9. shares_xts<-as.xts(shares[2:6],order.by=shares$time)11. times<-shares$time10.drawing<-function11. times<-shares$time12. if(as.Date(date)%in%times){13. 13. 15. share_date<-as.data.frame(share_date)15. share_date<-as.data.frame(share_date)16. 17. 17. share<-as.xts(share_date[2:6],order.by=share_date$time)19. dn.col="green",theme=chartTheme("white"))18. 19. dn.col="green",theme=chartTheme("white"))20. addMACD()21. 21. }else{print("Error")}22.}23.drawing("2021/8/12")#调用函数model_dataread.csv("G:/linear_model.csv"24.model_dataread.csv("G:/linear_model.csv"27.attach(model_data)28. 28. 30.x8<-0.395*x1+0.328*x2+0.426*x3-0.385*x4-3.433*x5-0.351*x629.30.x8<-0.395*x1+0.328*x2+0.426*x3-0.385*x4-3.433*x5-0.351*x631. 32.32.x7<-(x7-mean(x7))/sd(x7)34. 34. 35.y2<-y136.36.y2[y1<0]<-038.detach(model_data)37.y2[y1>0]<-38.detach(model_data)39.40.40.#建模42.model_data$x7<-x741. 42.model_data$x7<-x743.model_data$x8<-x844.44.model_data$y2<-y246.45.table(model_data$y2)#y值46.48.model1<-glm(y2~x1+x2+x3+x4+x5+x6,47. 48.model1<-glm(y2~x1+x2+x3+x4+x5+x6,49. data=model_data,family=binomial(link="logit"))50.50.summary(model1)52.51.AIC(model1);BIC(model1)52.53.model2<-glm(y2~x7+x8,data=model_data,54. 54. family=binomial(link="logit"))56.AIC(model2);BIC(model2)55.56.AIC(model2);BIC(model2)57.58.58.library(glmnet)60. family="binomial")5