精华内容
下载资源
问答
  • R中的统计模型

    千次阅读 2014-04-10 18:07:14
    R中的统计模型  这一部分假定读者已经对统计方法,特别是回归分析和方差分析有一定的了解。后面我们还会假定读者对广义线性模型和非线性模型也有所了解。R已经很好地定义了统计模型拟合中的一些前提条件,因此...

    R中的统计模型

     这一部分假定读者已经对统计方法,特别是回归分析和方差分析有一定的了解。后面我们还会假定读者对广义线性模型和非线性模型也有所了解。R已经很好地定义了统计模型拟合中的一些前提条件,因此我们能构建出一些通用的方法以用于各种问题。R提供了一系列紧密联系的统计模型拟合的工具,使得拟合工作变得简单。正如我们在绪论中提到的一样,基本的屏幕输出是简洁的,因此用户需要调用一些辅助函数来提取细节的结果信息。

     

    1定义统计模型的公式

    下面统计模型的模板是一个基于独立的方差齐性数据的线性模型

     用矩阵术语表示,它可以写成

           

    其中y是响应向量,X是模型矩阵(model matrix)或者设计矩阵(design ma-

    trix)。X的列 是决定变量(determiningvariable)。通常,列都是1,用来定义截距(intercept)项。

     例子

    在给予正式的定义前,举一些的例子可能更容易了解全貌。

    假定y,x,x0,x1,x2,...是数值变量,X是一个矩阵,而A,B,C,...是因子。下

    面的例子中,左边给出公式,右边给出该公式的统计模型的描述。

    y~x

    y~1+x

    二者都反映了y对x的简单线性模型。第一个公式包含了一个隐式的截距项,而第二个则是一个显式的截距项。

    y~0+x

    y~-1+x

    y~x-1                y对x过原点的简单线性模型(也就是说,没有截距项)。

    log(y)~x1+x2          y的变换形式log(y)对x1和x2进行的多重回归(有一个隐式的截距项)。

    y~poly(x,2)

    y~1+x+I(x^2)         y对x的二次多项式回归。第一种是正交多项式(orthogonal polynomial),第二种则显式地注明各项的幂次。

    y~X+poly(x,2)          y利用模型矩阵X和二次多项式项x进行多重回归。

    y~A                 y的单因素方差分析模型,类别由A决定。

    y~A+x               y的单因素协方差分析模型,类别由A决定,协方差项为x。

    y~A*B

    y~A+B+A:B

    y~B%in%A

    y~A/B                y对A和B的非可加两因子方差分析模型(two factor non-additive model)。前两个公式表示相同的交叉分类设计(crossedclassification),后两个公式表示相同的嵌套分类设计(nested classification)。抽象一点说,这四个公式指明同一个模型子空间。

    y~(A+B+C)^2

    y~A*B*C-A:B:C        三因子实验。该模型包括一个主效应(main effects)和两个因子的交互效应(interactions)。这两个公式等价。

    y~A*x

    y~A/x

    y~A/(1+x)-1            在A的各个水平独立拟合y对x的简单线性回归。三个公式的编码不一样。最后一个公式会对A各个水平分别估计截距项和斜率项的。

    y~A*B+Error(C)        一个实验设计有两个处理因素A和B以及因子C决定的误差分层(errorstrata)。如在裂区实验设计(split plotexperiment)中,所有区组(还包括子区组)都由因子C决定的。

    操作符~用来定义R的模型公式(model formula)。一个普通的线性模型式可以表示为

    response~op 1 term 1 op 2 term 2 op 3 term3...

    其中response是一个作为响应变量的向量或者矩阵,或者是一个值为向量/矩阵的表达式。op i是一个操作符。它要么是+要么是-,分别表示在一个模型中加入或者去掉某一项(公式第一项的操作符可选)。term i可以(1)是一个向量,矩阵表达式或者1,(2)因子,(3)是一个由因子,向量或矩阵通过公式操作符连接产生的公式表达式(formula expression)。

    基本上,公式中的项决定了模型矩阵中的列要么被加入要么被去除。1表示截距项,并且默认就已加入模型矩阵,除非显式地去除这一选项。

    公式操作符(formula operators)在效果上和用于程序Glim和Genstat中的Wilkinson&Rogers标记符(notation)相似。一个不可避免的改变是操作符.在R里面变成了:,因为点号在R里面是合法的命名字符。

    这些符号总结如下(参考Chambers&Hastie,1992,p.29):

    Y~M                       Y由模型M解释。

    M 1+M2                    同时包括M 1和M 2项。

    M 1-M2                    包括M 1但排除M 2项。

    M 1:M2                    M1和M 2的张量积(tensor product)。如果两项都是因子,那么将产生“子类”因子(subclasses factor)。

    M 1%in%M2                和M 1:M 2类似,但编码方式不一样。

    M 1*M 2                   M1+M 2+M 1:M 2.

    M 1/M 2                    M1+M 2%in%M1.

    M^n                        M的所有各项以及所有到n阶为止的“交互作用”项

    I(M)                       隔离M。M内所有操作符当一般的运算符处理。并且该项出现在模型矩阵中。

    注意,在常常用来封装函数参数的括弧中的操作符按普通的四则运算法则解

    释。I()是一个恒等函数(identity function),它使得常规的算术运算符可以用在模型公式中。还要特别注意模型公式仅仅指定了模型矩阵的列项,暗含了对参数项的指定。在某些情况下可能不是这样,如非线性模型的参数指定。

     

    1对照

    我们至少要知道模型公式是如何指定模型矩阵的列项的。对于连续变量这是比较简单的,因为每一个变量对应于模型矩阵的一个列(如果模型中包含截距,会在矩阵中列出值都是1的一列)。

    对于一个k-水平的因子A该如何处理呢?无序和有序因子给出的结论是不一样的。对于无序因子,因子第2,...,第k不同水平的指标产生k?1列。(因此隐含的参数设置就是把其他水平和第一个水平的响应程度进行比较)。对于有序因子,k-1列是在1,...,k上的正交项(orthogonal polynomial),并且忽略常数项。

    尽管这里的回答有点复杂,但这不是事情的全部。首先在含有一个因子项的模型中忽略截距项,这一项将会被编入指示所有因子水平的k列中。其次整个行为可以通过options设置参数contrasts而改变。R的默认设置为

    options(contrasts=c("contr.treatment","contr.poly"))

    提这些内容的主要原因是R和S对无序因子采用不同的默认值。S采用Helmert对照。因此,当你需要比较你的结果和某本书上或论文上用SPLUS代码的结果时,你必须设置

    options(contrasts=c("contr.helmert","contr.poly"))

    这是一个经过认真考虑的改变。因为处理对照(treatment contrast)(R默认)对于新手是比较容易理解的。

    这还没有结束,因为在各个模型的各个项中对照方式可以用函数contrasts和C重新设置。

    我们还没有考虑交互作用项:这些交互作用项将会产生各分量项的乘积。

    尽管细节是复杂的,R里面的模型公式在要求不是太离谱的情况下可以产生统计专家所期望的各种模型。提供模型公式的各种扩展特性是让R更灵活。例如,利用关联项而非主要效应的模型拟合常常会产生令人惊讶的结果,不过这些仅仅为统计专家们设计的。

     

    2线性模型

    对于常规的多重模型(multiple model)拟合,最基本的函数是lm()。下面是调

    用它的方式的一种改进版:

    >fitted.model<-lm(formula,data=data.frame)

    例如

    >fm2<-lm(y~x1+x2,data=production)

    将会拟合y对x1和x2的多重回归模型(和一个隐式的截距项)。

    一个重要的(技术上可选)参数是data=production。它指定任何构建这个模型的变量首先必须来自数据框production。这里不需要考虑数据框production是否被绑定在搜索路径中。

     

    3提取模型信息的泛型函数

    lm()的返回值是一个模型拟合结果对象;技术上就是属于类"lm"的一个结果列表。关于拟合模型的信息可以用适合对象类"lm"的泛型函数显示,提取,图示等等。这包括

    add1 coef effects kappa predict residuals

    alias deviance family labels print step

    anova drop1 formula plot proj summary

    其中一些常用的泛型函数可以简洁描述如下。

    anova(object 1,object2)     比较一个子模型和外部模型,并且产生方差分析表。

    coef(object)               提取回归系数(矩阵)。全称:coefficients(object).

    deviance(object)           残差平方和,若有权重可加权。

    formula(object)            提取模型公式信息。

    plot(object)               产生四个图,显式残差,拟合值和一些诊断图。

    predict(object,newdata=data.frame)提供的数据框必须有同原始变量一样标签的变量。结果是对应于data.frame中决定变量预测值的向量或矩阵。

    predict.gam(object,

    newdata=data.frame)        predict.gam()是安全模式的predict()。它可以用于lm,glm和gam拟合对象。在正交多项式作为原始的基本函数并且增加新数据意味着必须使用不同的原始基本函数。

    print(object)                简要打印一个对象的内容。常常隐式使用。

    residuals(object)            提取残差(矩阵),有权重时可加权,省略方式:resid(object)。

    step(object)                 通过增加或者减少模型中的项并且保留层次来选择合适的模型。在逐步搜索过程中,AIC(Akaike信息规范)值最大的模型将会被返回。

    summary(object)             显示较详细的模型拟合结果。

     

    4方差分析和模型比较

    aov(formula,data=data.frame)和函数lm()非常的相似,在泛型函数提取模型信息部分列出的泛型函数同样适用。

    需要注意的是aov()还允许分析多误差层次(multiple error strata)的模型,如

    裂区实验设计(split plot experiments),利用区组内信息进行的平衡不完全区组设

    计(balanced incomplete block design)等。模型公式

    response mean.formula+Error(strata.formula)

    指定了一个多层次实验设计,误差层由strata.formula定义。最简单的情况是,strata.formula是单因素的。它定义了一个双层次的实验,也就是研究在这些因子的水平内或者水平间的实验响应。

    例如,已知所有的决定变量因子,模型公式可以设计如下:

    >fm<-aov(yield~v+n*p*k+Error(farms/blocks),data=farm.data)

    这常常用来描述一个同时含有均值模型v+n*p*k和三个误差层次(“农田之间”,“农田内但在区组之间”和“区组内”)的实验。

     

    方差表的分析实际上是为拟合模型序列(sequence)进行的。在模型序列的特定地方增加特定的项会使残差平方和降低。因此仅仅在正交实验中,模型中增加项的次序是没有影响的。

    在多层实验设计(multistratum experiments)中,程序首先把响应值依次投射到各个误差层次上,并且用均值模型去拟合各个投射。细节内容可以参考Chambers&Hastie(1992)。

    除了使用常规的方差分析表(ANOVA table),你还可以直接用函数anova()来比较两个模型。这种方法更为灵活。

    >anova(fitted.model.1,fitted.model.2,...)

    结果将是一个方差分析表以显示依次加入的拟合模型的差异。需要比较的拟合模

    型常常是等级序列(hierarchical sequence)。这个和默认的设置实际上没有差别,只是使它更容易理解和控制。

     

    5更新拟合模型

    函数update()是一个非常便利的函数。它允许拟合一个比原先模型增加或减少一个项的模型。它的形式是

    >new.model<-update(old.model,new.formula)

    在new.formula中,公式中包含的句点,.,仅仅表示“旧的公式模型中的对应部

    分”。例如

    >fm05<-lm(y~x1+x2+x3+x4+x5,data=production)

    >fm6<-update(fm05,.~.+x6)

    >smf6<-update(fm6,sqrt(.)~.)

    这将分别拟合从数据框production中得到的五个变量的多重回归模型,拟合额外增加一个变量的六个回归量的模型,和进一步对响应值进行平方根变换后的模型拟合。

    注意参数data=在最开始调用模型拟合函数的时候指定。这个信息将会通过拟合模型对象传递给函数update()及其相关者。

    符号.同样可以用在其他情况下,不过含义有点不同。如

    >fmfull<-lm(y~.,data=production)

    它将会拟合响应量y对数据框production中所有变量回归的模型。

    其他研究逐步回归的函数是add1(),drop1()和step()。从字面上就可以看出这些函数的意义,更细节的内容可以参考在线帮助文档。

     

    6广义线性模型

    广义线性建模是线性建模的一种发展,它通过一种简洁而又直接的方式使得线性模型既适合非正态分布的响应值又可以进行线性变换。广义线性模型是基于下面一系列假设前提的:

    (1)有一个有意思的响应变量y和一系列刺激变量(stimulusvariable)x1,x2,...。

    这些刺激变量决定响应变量的最终分布。

    (2)刺激变量仅仅通过一个线性函数影响响应值y的分布。该线性函数称为线性预测器(linear predictor),常常写成

    η=β1x12x2+···+βpxp,

    因此xi当且仅当βi=0时对y的分布没有影响。

    (3)y分布的形式为

     

     

    其中是尺度参数(scale parameter)(可能已知),对所有观测恒定,A是一个先验的权重,假定知道但可能随观测不同有所不同,μ是y的均值。也就是说假定y的分布是由均值和一个可能的尺度参数决定的。

    (4)均值μ是线性预测器的平滑可逆函数(smooth invertible function):

    μ=m(η),η=m-1(μ)=l(μ)

    其中的反函数(inverse function)l()被称为关联函数(link function)。

    这些假定比较宽松,足以包括统计实践中大多数有用的统计模型,同时也足够严谨,使得可以发展参数估计和统计推论(estimation and inference)中一致的方法(至少可以近似一致)。读者如果想了解这方面最新的进展,可以参考McCullagh&Nelder(1989)或者Dobson(1990)。

    6.1族

    R提供了一系列广义线性建模工具,从类型上来说包括高斯(gaussian),二项式,泊松(poisson),逆高斯(inverse gaussian)和伽马(gamma)模型的响应变量分布以及响应变量分布无须明确给定的拟似然(quasi-likelihood)模型。在后者,方差函数(variance function)可以由均值的函数指定,但在其它情况下,该函数可以由响应变量的分布得到。每一种响应分布允许各种关联函数将均值和线性预测器关联起来。这些自动可用的关联函数如下表所示:

    族名字                      关联函数

    Binomial                logit,probit,log,cloglog

    Gaussian                identity,log,inverse

    Gamma                 identity,inverse,log

    inverse.aussian          1/mu^2,identity,inverse,log

    poisson                identity,log,sqrt

    quasi                  logit,probit,cloglog,identity,

    inverse,log,1/mu^2,sqrt

    这些用于模型构建过程中的响应分布,关联函数和各种其他必要的信息统称为广义线性模型的族(family)。

    6.2 glm()函数

    既然响应的分布仅仅通过单一的一个线性函数依赖于刺激变量,那么用于线性模型的机制同样可以用于指定一个广义模型的线性部分。但是族必须以一种不同的方式指定。

    R用于广义线性回归的函数是glm(),它的使用形式为

    >fitted.model<-glm(formula,family=family.generator,data=data.frame)

    和lm()相比,唯一的一个新特性就是描述族的参数family.generator。它其实是一个函数的名字,这个函数将产生一个函数和表达式列表用于定义和控制模型的构建与估计过程。尽管这些内容开始看起来有点复杂,但它们非常容易使用。

    这些名字是标准的。程序给定的族生成器可以参见族部分表格中的“族名”。当选择一个关联函数时,该关联函数名和族名可以同时在括弧里面作为参数设定。在拟(quasi)家族里面,方差函数也是以这种方式设定。

    一些例子可能会使这个过程更清楚。

    gaussian族

    命令

    >fm<-glm(y~x1+x2,family=gaussian,data=sales)

    和下面的命令结果一致

    >fm<-lm(y~x1+x2,data=sales)

    但是效率上,前者差一点。注意,高斯族没有自动提供关联函数设定的选项,因此不允许设置参数。如一个问题需要用非标准关联函数的高斯族,那么只能采用我们后面讨论的拟族。

    二项式族

    考虑Silvey(1970)提供的一个人造的小例子。

    在Kalythos的Aegean岛上,男性居民常常患有一种先天的眼科疾病,并且随着年龄的增长而变的愈明显。现在搜集了各种年龄段岛上男性居民的样本,同时记录了盲眼的数目。数据展示如下:

    Age:20 35 45 55 70

    No.:tested:50 50 50 50 50

    No.:blind:6 17 26 37 44

    我们想知道的是这些数据是否吻合logistic和probit模型,并且分别估计各个模型的LD50,也就是一个男性居民盲眼的概率为50%时候的年龄。

    如果y和n分别是年龄为x时的盲眼数目和检测样本数目,两种模型的形式都为

    y~B(n,F(β01x))

    其中在probit模型中,F(z)=Φ(z)是标准的正态分布函数,而在logit模型(默认)中,F(z)=ez/(1+ez)。这两种模型中,

    LD50=-β0/β1

    ,即分布函数的参数为0时所在的点。

    第一步是把数据转换成数据框,

    >kalythos<-data.frame(x=c(20,35,45,55,70),n=rep(50,5),y=c(6,17,26,37,44))

    在glm()拟合二项式模型时,响应变量有三种可能性:

    (1)如果响应变量是向量,则假定操作二元(binary)数据,因此要求是0/1向量。

    (2)如果响应变量是双列矩阵,则假定第一列为试验成功的次数第二列为试验失败

    的次数。

    (3)如果响应变量是因子,则第一水平作为失败(0)考虑而其他的作为‘成功’(1)考虑。

    这里,我们采用的是第二种惯例。我们在数据框中增加了一个矩阵:

    >kalythos$Ymat<-cbind(kalythos$y,kalythos$n-kalythos$y)

    为了拟合这些模型,我们采用

    >fmp<-glm(Ymat~x,family=binomial(link=probit),data=kalythos)

    >fml<-glm(Ymat~x,family=binomial,data=kalythos)

    既然logit的关联函数是默认的,因此我们可以在第二条命令中省看拟合结果,我们使用

    >summary(fmp)

    >summary(fml)

    两种模型都拟合的很好。为了计算LD50,我们可以利用一个简单

    >ld50<-function(b)-b[1]/b[2]

    >ldp<-ld50(coef(fmp));ldl<-ld50(coef(fml));c(ldp,ldl)

    从这些数据中得到的年龄分别是43.663年和43.601年。

    Poisson模型

    Poisson族默认的关联函数是log。在实际操作中,这一族常常用于拟合计数资料的Poisson对数线性模型。这些计数资料的实际分布往往符合二项式分布。这是一个非常重要而又庞大的话题,我们不想在这里深入展开。它甚至是非-高斯广义模型内容的主要部分。

    有时候,实践中产生的Poisson数据在对数或者平方根转化后可当作正态数据处理。作为后者的另一种选择是,一个Poisson广义线性模型可以通过下面的方式拟合:

    >fmod<-glm(y~A+B+x,family=poisson(link=sqrt),data=worm.counts)

    拟似然模型

    对于所有的族,响应变量的方差依赖于均值并且拥有作为乘数(multiplier)的尺度参数。方差对均值的依赖方式是响应分布的一个特性;例如对于poisson分布Var[y]=μ。

    对于拟似然估计和推断,我们不是设定精确的响应分布而是设定关联函数和方差函数的形式,因为关联函数和方差函数都依赖于均值。既然拟似然估计和gaussian分布使用的技术非常相似,因此这一族顺带提供了一种用非标准关联函数或者方差函数拟合gaussian模型的方法。

    例如,考虑非线性回归的拟合

    y=θ1z1/(z22)+e

    同样还可以写成

    y=1/(β1x12x2)+e

    其中x1=z2/z1,x2=-1/x11=1/θ1 andβ221。假如有适合的数据框,我们可以如下进行非线性拟合

    >nlfit<-glm(y~x1+x2-1,

    family=quasi(link=inverse,variance=constant),

    data=biochem)

    如果需要的话,读者可以从其他手册或者帮助文档中得到更多的信息。

     

    7非线性最小二乘法和最大似然法模型

    特定形式的非线性模型可以通过广义线性模型(glm())拟合。但是许多时候,我们必须把非线性拟合的问题作为一个非线性优化的问题解决。R的非线性优化程序是optim(),nlm()和nlminb()(自R2.2.0开始)。二者分别替换SPLUS的ms()和nlminb()但功能更强。我们通过搜寻参数值使得缺乏度(lack-of-fit)指标最低,如nlm()就是通过循环调试各种参数值得到最优值。和线性回归不同,程序不一定会收敛到一个稳定值。nlm()需要设定参数搜索的初始值,而参数估计是否收敛在很大程度上依赖于初始值设置的质量。

    7.1最小二乘法

    拟合非线性模型的一种办法就是使误差平方和(SSE)或残差平方和最小。如果观测到的误差极似正态分布,这种方法是非常有效的。

    下面是例子来自Bates&Watts(1988),51页。具体数据是:

    >x<-c(0.02,0.02,0.06,0.06,0.11,0.11,0.22,0.22,0.56,0.56,1.10,1.10)

    >y<-c(76,47,97,107,123,139,159,152,191,201,207,200)

    被拟合的模型是:

    >fn<-function(p)sum((y-(p[1]*x)/(p[2]+x))^2)

    为了进行拟合,我们需要估计参数初始值。一种寻找合理初始值的办法把数据图形化,然后估计一些参数值,并且利用这些值初步添加模型曲线。

    >plot(x,y)

    >xfit<-seq(.02,1.1,.05)

    >yfit<-200*xfit/(0.1+xfit)

    >lines(spline(xfit,yfit))

    当然,我们可以做的更好,但是初始值200和0.1应该足够了。现在做拟合:

    >out<-nlm(fn,p=c(200,0.1),hessian=TRUE)

    拟合后,out$minimum是误差的平方和(SSE),out$estimate是参数的最小二乘估计值。为了得到参数估计过程中近似的标准误(SE),我们可以:

    >sqrt(diag(2*out$minimum/(length(y)-2)*solve(out$hessian)))

    上述命令中的2表示参数的个数。一个95%的信度区间可以通过±1.96 SE计算得到。我们可以把这个最小二乘拟合曲线画在一个新的图上:

    >plot(x,y)

    >xfit<-seq(.02,1.1,.05)

    >yfit<-212.68384222*xfit/(0.06412146+xfit)

    >lines(spline(xfit,yfit))

    标准包stats提供了许多用最小二乘法拟合非线性模型的扩充工具。我们刚刚拟合过的模型是Michaelis-Menten模型,因此可以利用下面的命令得到类似的结论。

    >df<-data.frame(x=x,y=y)

    >fit<-nls(y~SSmicmen(x,Vm,K),df)

    >fit

    Nonlinear regression model

    model:         y~SSmicmen(x,Vm,K)

    data:          df

    Vm                K

    212.68370711      0.06412123

    residual sum-of-squares:    1195.449

    >summary(fit)

    Formula:  y~SSmicmen(x,Vm,K)

    Parameters:

    Estimate   Std.Error    tvalue      Pr(>|t|)

    Vm      2.127e+02 6.947e+00   30.615     3.24e-11

    K       6.412e-02  8.281e-03    7.743     1.57e-05

    Residual standard error: 10.93 on 10degrees of freedom

    Correlation of Parameter Estimates:

    Vm

    K 0.7651

    7.2最大似然法

    最大似然法(Maximum likelihood)也是一种非线性拟合方法。它甚至可以用在误差非正态的数据中。这种方法估计的参数将会使得对数似然值最大或者负的对数似然值最小。下面的例子来自Dobson(1990),pp.:108–111。这个例子对剂量-响应数据拟合logistic模型(当然也可以用glm()拟合)。数据是:

    >x<-c(1.6907,1.7242,1.7552,1.7842,1.8113,1.8369,1.8610,1.8839)

    >y<-c(6,13,18,28,52,53,61,60)

    >n<-c(59,60,62,56,63,59,62,60)

    要使负对数似然值最小,则:

    >fn<-function(p)

    sum(-(y*(p[1]+p[2]*x)-n*log(1+exp(p[1]+p[2]*x))+log(choose(n,y))))

    我们选择一个适当的初始值,开始拟合:

    >out<-nlm(fn,p=c(-50,20),hessian=TRUE)

    拟合后,out$minimum就是负对数似然值,out$estimate就是最大似然拟合的参数值。为了得到拟合过程近似的标准误,我们可以:

    >sqrt(diag(solve(out$hessian)))

    参数估计的95%信度期间可由估计值±1.96 SE计算得到。

     

    8一些非标准模型

    我们简单提一下R里面某些用于某些特殊回归和数据分析问题的工具。

    (1)混合模型(Mixed models)。用户捐献包nlme里面提供了函数lme()和nlme()。这些函数可以用于混合效应模型(mixed-effects models)的线性和非线性回归。也就是说在线性和非线性回归中,一些系数受随机因素的影响。这些函数

    需要充分利用公式来描述模型。

    (2)局部近似回归(Localapproximating regressions)。函数loess()利用局部加权回归进行一个非参数回归。这种回归对显示一组凌乱数据的趋势和描述大数据集的整体情况非常有用。函数loess和投影跟踪回归(projection pursuit regression)的代码一起放在标准包stats中。

    (3)稳健回归(Robustregression)。有多个函数可以用于拟合回归模型,同时尽量不受数据中极端值的影响。在推荐包MASS中的函数lqs为高稳健性的拟合提供了最新的算法。另外,稳健性较低但统计学上高效的方法同样可以在包MASS中得到,如函数rlm。

    (4)累加模型(Additive models)。这种技术期望可以通过决定变量的平滑累加函数(smooth additive function)构建回归函数。一般来说,每个决定变量都有一个平滑累加函数。用户捐献的包acepack里面的函数avas和ace以及包mda里面的函数bruto和mars为这种技术提供了一些例子。这种技术的一个扩充是用户捐献包gam和mgcv里面实现的广义累加模型。

    (5)树型模型(Tree-basedmodels)。除了利用外在的全局线性模型预测和解释数据,还可以利用树型模型递归地在决定性变量的判断点上将数据的分叉分开。这样做会把数据最终分成几个不同组,使得组内尽可能相似而组间尽可能差异。这样常常会得到一些其他数据分析方法不能产生的的信息。模型可以用一般的线性模型形式指定。该模型拟合函数是tree(),而且许多泛型函数,如plot()和text()都可以很好的用于树型模型拟合结果的图形显示。R里面的树型模型函数可以通过用户捐献的包rpart和tree得到。

     

    展开全文
  • R中GARCH模型

    2017-12-28 23:31:45
    R中建立GARCH(1,1)模型的代码,帮助大家解决时间序列问题
  • 风险评估模型蒙特卡洛模型Note from Towards Data Science’s editors: While we allow independent authors to publish articles in accordance with our rules and guidelines, we do not endorse each author’s ...

    风险评估模型蒙特卡洛模型

    Note from Towards Data Science’s editors: While we allow independent authors to publish articles in accordance with our rules and guidelines, we do not endorse each author’s contribution. You should not rely on an author’s works without seeking professional advice. See our Reader Terms for details.

    Towards Data Science编辑的注意事项: 尽管我们允许独立作者按照我们的 规则和指南 发表文章 ,但我们不认可每位作者的贡献。 您不应在未征求专业意见的情况下依赖作者的作品。 有关 详细信息, 请参见我们的 阅读器条款

    A common problem when evaluating a portfolio manager is that the history of returns is often so short that estimates of risk and performance measures can be highly unreliable. A similar problem occurs when testing a new trading strategy. Even if you have a fairly long history for the strategy’s performance, often you only have observations over a single market cycle which can make it difficult to evaluate how your strategy would have held up in other markets. If you trade stocks you have probably heard the refrain: “I’ve never seen a bad back test”.

    在评估投资组合经理时,一个常见的问题是收益的历史通常很短,以至于风险和绩效指标的估计可能非常不可靠。 测试新的交易策略时会发生类似的问题。 即使您对策略的执行历史已有很长的历史,但通常您只能在单个市场周期中观察到,这可能使您难以评估您的策略在其他市场中的表现。 如果您交易股票,您可能会听到这样的话:“我从未见过糟糕的回测”。

    One method to address this deficiency is through Factor Model Monte Carlo (FMMC). FMMC can be used to estimate a factor model based on a set of financial and economic factors that reliably explain the returns of the fund manager. We can then simulate returns to determine how the manager would have performed in a wide variety of market environments. The end result is a model that produces considerably better estimates for risk and performance than if we simply used the return series available to us.

    解决此缺陷的一种方法是通过因素模型蒙特卡洛(FMMC) 。 FMMC可用于基于一组可靠地解释基金经理收益的财务和经济因素来估计因素模型。 然后,我们可以模拟收益,以确定经理在各种市场环境中的表现。 最终结果是,与仅使用可用的收益序列相比,该模型可以更好地估算风险和绩效。

    任务和设置 (The Task and Set Up)

    For this case study, we will be analyzing the returns for the new hedge fund Aric’s Hedge Fund; hereafter known as AHF. The hedge fund case is particularly interesting because hedge funds can use leverage, invest in any asset class, go long or short, and use many different instruments. Hedge funds are often very secretive about their strategy and holdings. Thus, having a reliable risk model to explain the source of their returns is essential.

    在本案例研究中,我们将分析新对冲基金Aric's Hedge Fund的回报; 以下简称AHF 。 对冲基金的案例特别有趣,因为对冲基金可以利用杠杆,投资于任何资产类别,做多或做空以及使用许多不同的工具。 对冲基金通常对其策略和持股非常保密。 因此,拥有可靠的风险模型来解释其收益的来源至关重要。

    Keep in mind that Aric’s Hedge Fund is not a real hedge fund (I’m Aric, I don’t have a hedge fund), but this is a real series of returns. I obtained the returns for a hedge fund in operation that we invest in where I work so the results of this study are applicable to a real-world scenario.

    请记住,Aric的对冲基金不是真正的对冲基金(我是Aric,我没有对冲基金), 但这是一系列真实的回报。 我获得了对冲基金的回报,该对冲基金投资于我工作的地方,因此本研究的结果适用于现实情况。

    We have data for Aric’s Hedge Fund from January 2010 to March 2020. For the purpose of this post and evaluating the accuracy of our model we will pretend as though AHF is pretty new to the scene and that we only have data from January 2017 through March 2020. To overcome the data deficiency, we will build a factor model on the basis of this “observed data” and then utilize the entire data series to evaluate the accuracy of our simulation for assessing the risk and performance statistics.

    我们拥有2010年1月至2020年3月Aric对冲基金的数据。出于这篇文章的目的,并评估我们模型的准确性,我们将假装 AHF尚不成熟,并且我们仅提供2017年1月至3月的数据。 2020年。为克服数据不足,我们将在此“观察数据”的基础上建立一个因子模型,然后利用整个数据系列来评估模拟的准确性,以评估风险和绩效统计数据。

    The below graph shows the cumulative return of AHF since January 2010. The data to the right of the red line represents the “observed period”.

    下图显示了自2010年1月以来AHF的累计收益。红线右侧的数据代表“观察期”。

    Image for post

    We will be conducting the analysis in R using the extensive library of packages available therein including: PerformanceAnalytics and quantmod. Aside from the hedge fund returns series, all of the factor data can be obtained freely from Yahoo! Finance, the Federal Reserve Bank of St. Louis FRED Database, and the Credit Suisse Hedge Fund Indices; you have to sign up for Credit Suisse to access the indices, but still…free.

    我们将使用其中提供的广泛的软件包库在R中进行分析,包括: PerformanceAnalyticsquantmod 。 除了对冲基金收益系列以外,所有因子数据都可以从Yahoo!免费获得。 金融 ,圣路易斯联邦储备银行FRED数据库瑞士信贷对冲基金指数 ; 您必须注册瑞士信贷才能访问该指数,但仍然…免费。

    模型估计 (Model Estimation)

    A common technique in empirical finance is to explain changes in asset prices based on a set of common risk factors. The simplest and most well-known factor model is the Capital Asset Pricing Model (CAPM) of William Sharpe. The CAPM is specified as follows:

    经验金融中的一种常见技术是基于一组常见风险因素来解释资产价格的变化。 最简单和最著名的因素模型是William Sharpe的资本资产定价模型(CAPM)。 CAPM指定如下:

    Image for post

    Where:

    哪里:

    • ri = Return of asset ‘i’

      ri =资产“ i”的收益
    • m = Return of market index ‘m’

      m =市场指数“ m”的回报
    • ∝= Excess return

      ∝ =超额收益
    • ẞ = Exposure to the Market Risk factor

      ẞ=暴露于市场风险因素
    • Ɛi = Idiosyncratic error term

      = i =特异误差项

    Market risk, or “systematic” risk, serves as a summary measure for all of the risks to which financial assets are exposed. This may include recessions, inflation, changes in interest rates, political turmoil, natural disasters, etc. Market risk is usually proxied by the returns on a large index like the S&P 500 and cannot be reduced through diversification.

    市场风险或“系统性”风险,是金融资产​​所面临的所有风险的汇总度量。 这可能包括衰退,通货膨胀,利率变化,政治动荡,自然灾害等。市场风险通常由标普500指数等大型指数的回报所替代,并且无法通过多元化降低。

    ẞ (i.e. Beta) represents an asset’s exposure to market risk. A Beta = 1 would imply that the asset is as risky as the market, Beta >1 would imply more risk than the market, while a Beta < 1 would imply less risk.

    ẞ(即Beta)代表资产的市场风险。 Beta = 1表示资产的风险与市场一样,Beta> 1的风险比市场高,而Beta <1的风险比市场低。

    Ɛ is idiosyncratic risk and represents the portion of the return that cannot be explained by the Market Risk factor.

    Ɛ是特质风险,代表不能由市场风险因素解释的收益部分。

    We will extend the CAPM to include additional risk factors which the literature have shown to be important for explaining asset returns. Aric’s Hedge Fund runs a complicated strategy using many different asset classes and instruments so it’s certainly plausible that it would be exposed to a broader set of risks beyond the traditional market index. The general form of our factor model is as follows:

    我们将把CAPM扩展到包括其他风险因素,这些因素已被文献证明对于解释资产收益很重要。 Aric的对冲基金使用许多不同的资产类别和工具来执行一项复杂的策略,因此,它可能会面临比传统市场指数更广泛的风险。 我们的因子模型的一般形式如下:

    Image for post

    All the above says is that returns (r) are explained by a set of risk factors j=1…k where r j is the return for factor ‘j’ and ẞ j is the exposure. Ɛ is the idiosyncratic error. Thus, if we can estimate the ẞ j, then we can leverage the long history of factor returns (r j) calculate conditional returns for AHF. Finally, if we can reasonably estimate the distribution of Ɛ then we can build randomness into AHF’s return series. This enables us to fully capture the variety of returns that we could observe.

    上面所说的全部是,回报(r)由一组风险因子j = 1…k解释,其中rj是因子“ j”的回报,而ẞj是风险敞口。 Ɛ是特质错误。 因此,如果我们可以估计ẞj,那么我们就可以利用长期的要素收益率(rj)计算AHF的条件收益率。 最后,如果我们可以合理估计distribution的分布,则可以将随机性构建到AHF的收益序列中。 这使我们能够充分捕捉我们可以观察到的各种回报。

    The FMMC method will take place in three parts:

    FMMC方法将分为三个部分:

    • Part A: Data Acquisition, Clean Up and Processing

      A部分:数据采集,清理和处理
    • Part B: Model Estimation

      B部分:模型估算
    • Part C: Monte Carlo Simulation

      C部分:蒙特卡洛模拟

    A部分:数据采集,清理和处理 (Part A: Data Acquisition, Clean Up and Processing)

    For the factor model I will be using a set of financial and economic variables aimed at measuring different sources of risk and return. Again, all the data used in this study are freely available from Yahoo! Finance, the FRED Database, and Credit Suisse.

    对于因子模型,我将使用一组金融和经济变量,旨在衡量风险和回报的不同来源。 同样,本研究中使用的所有数据均可从Yahoo!免费获得。 金融,FRED数据库和瑞士信贷。

    We’ll begin with the FRED data. Next to each variable I have placed the unique identifier that you can query from the database.

    我们将从FRED数据开始。 在每个变量旁边,我放置了可以从数据库查询的唯一标识符。

    FRED Variables:

    FRED变量:

    • 5-Year Inflation Expectations, 5-Years Forward. (T5YIFRM)

      未来5年通胀预期,未来5年。 (T5YIFRM)
    • Term Spread: 10-Year minus 3-month Treasury Yield Spread. (T10Y3M)

      期限利差:10年期减去3个月的美国国债收益率利差。 (T10Y3M)
    • Credit spread premium: Moody’s Baa corp bond yield minus 10-year Treasury yield. (BAA10Y)

      信用利差溢价:穆迪的Baa公司债券收益率减去10年期美国国债收益率。 (BAA10Y)
    • 3-month T-bill rate. (DGS3MO)

      3个月期国库券利率。 (DGS3MO)
    • TED Spread: 3-Month LIBOR Minus 3-Month Treasury Yield. (TEDRATE)

      TED利差:3个月伦敦银行同业拆借利率减去3个月国库券收益率。 (TEDRATE)
    • International bond yield: 10-Year government bond yields for Euro Area. (IRLTLT01EZM156N)

      国际债券收益率:欧元区10年期政府债券收益率。 (IRLTLT01EZM156N)
    • Corporate Bond Total Return Index: ICE BofAML Corp bond master total return index; in levels. (BAMLCC0A0CMTRIV)

      公司债券总收益指数:ICE BofAML Corp债券主总收益指数; 在水平上。 (BAMLCC0A0CMTRIV)
    • CBOE Volatility Index (i.e. the VIX). (VIXCLS)

      CBOE波动率指数(即VIX)。 (VIXCLS)
    • CBOE Volatility Index of US 10-Year Treasuries (i.e. the Treasury VIX). (VXTYN)

      美国10年期美国国债(即美国国债VIX)的CBOE波动率指数。 (VXTYN)

    The FRED API leaves something to be desired and does not allow you to pull data in a consistent way. The returns of AHF are monthly so our model will need to be estimated using monthly data. However, FRED retrieves data at the highest available frequency so daily data always comes in as daily. Furthermore, the data is retrieved from the beginning of the series, so you end up getting a lot of NAs. As such, we will need to do a little clean up before we proceed.

    FRED API有一些不足之处,并且不允许您以一致的方式提取数据。 AHF的回报是每月的,因此我们的模型需要使用每月的数据进行估算。 但是,FRED会以最高可用频率检索数据,因此每日数据始终以每日形式出现。 此外,数据是从本系列的开始检索的,因此您最终会获得很多NA。 因此,在继续之前,我们需要进行一些清理。

    The following segments of R code show loading the identifiers into variables and separate queries to FRED for the daily and monthly data. The daily data is cleaned and converted to a monthly frequency. I’ve tried to comment the code as much as possible so you can see what’s happening.

    R代码的以下各节显示将标识符加载到变量中,并分别向FRED查询每日和每月数据。 每日数据将被清理并转换为每月一次。 我尝试了尽可能多地注释代码,以便您可以看到发生了什么。

    The FRED data is good to go. The other set of variables that we will need are financial market indices. Growth, Value and Size indices feature prominently in asset pricing models such as the Fama-French 3-Factor Model and I take the same approach here. Returns from financial indices are obtained from the venerable Yahoo! Finance.

    FRED数据很好。 我们将需要的另一组变量是金融市场指数。 增长,价值和规模指数在诸如Fama-French 3-Factor Model之类的资产定价模型中具有突出的地位,在此我采用相同的方法。 财务指数的收益来自古老的Yahoo! 金融。

    Yahoo! Finance Variables:

    雅虎! 财务变量:

    • Value: Russell 1000 Value Index (^RLV)

      值:罗素1000价值指数(^ RLV)
    • Growth: Russell 1000 Growth Index (^RLG)

      成长:罗素1000成长指数(^ RLG)
    • Size: Russell 2000 Index (^RUT)

      大小:罗素2000指数(^ RUT)
    • Market: S&P 500 Index (^GSPC)

      市场:标普500指数(^ GSPC)
    • International: MSCI EAFE (EFA)

      国际:MSCI EAFE(EFA)
    • Bonds: Barclays Aggregate Bond Index (AGG)

      债券:巴克莱综合债券指数(AGG)

    Lastly, we’ll load in the hedge fund specific indices courtesy of Credit Suisse (CS). Obtaining the index requires a few extra steps as the data needs to manually downloaded to Excel from the Credit Suisse website and then loaded into R. Each index corresponds to a specific hedge fund strategy.

    最后,我们将根据瑞士信贷(CS)的数据加载对冲基金的特定指数。 获取索引需要一些额外的步骤,因为需要将数据从瑞士信贷网站手动下载到Excel,然后加载到R中。每个索引对应于特定的对冲基金策略。

    Credit Suisse Variables:

    瑞士信贷变量:

    • Convertible Bond Arbitrage Index (CV_ARB)

      可转换债券套利指数(CV_ARB)
    • Emerging and Frontier Markets Index (EM_MRKT)

      新兴和前沿市场指数(EM_MRKT)
    • Equity Market Neutral Index (EQ_Neutral)

      股市中性指数(EQ_Neutral)
    • Event Driven Index (EVT_DRV)

      事件驱动索引(EVT_DRV)
    • Distressed Opportunities Index (DISTRESS)

      苦恼机会指数(DISTRESS)
    • Multi-Strategy Event Driven Index (MS_EVT)

      多策略事件驱动索引(MS_EVT)
    • Event Driven Risk Arbitrage Index (RISK_ARB)

      事件驱动风险套利指数(RISK_ARB)
    • Fixed Income Arbitrage Index (FI_ARB)

      固定收益套利指数(FI_ARB)
    • Global Macro Index (GL_MACRO)

      全球宏观指数(GL_MACRO)
    • Equity Long-Short Index (EQ_LS)

      股票多空指数(EQ_LS)
    • Managed Futures Index (MNGD_FT)

      管理期货指数(MNGD_FT)
    • Multi-Strategy Hedge Fund Index (MS_HF)

      多策略对冲基金指数(MS_HF)

    B部分:模型估算 (Part B: Model Estimation)

    Recall that for the purpose of this case study we are “pretending” as though we only have data for AHF from January 2017 through March 2020 (i.e. the sample period). In reality we have data going back to January 2010. We will use the data in the sample period to calibrate the factor model and then compare the results from the simulation to the long-run risk and performance over the full period of January 2010 to March 2020.

    回想一下,就本案例研究而言,我们“假装”为好像我们只有从2017年1月到2020年3月(即采样期)的AHF数据。 实际上,我们拥有的数据可以追溯到2010年1月。我们将使用采样期间的数据来校准因子模型,然后将模拟结果与2010年1月至3月整个期间的长期风险和绩效进行比较2020年。

    Model estimation has 2-steps:

    模型估算分为两步:

    1. Estimate a Factor Model: Using the common “short” history of asset and factor returns, compute a factor model with intercept, factor betas j=1…k, and residuals

      估计因子模型:使用常见的资产和因子收益的“短”历史记录,计算具有截距,因子beta j = 1…k和残差的因子模型

    2. Estimate Error Density: Use the residuals from the factor model to fit a suitable density function from which we can draw.

      估计误差密度:使用因子模型中的残差来拟合一个合适的密度函数,我们可以从中得出该密度函数。

    I have proposed 27 risk factors to explain the returns of AHF, but I don’t know ahead of time which form the best prediction. It could be that some factors are irrelevant and reduce the explanatory power of the model. In order to select an optimal model, I use an Adjusted-R2 based best-subset selection algorithm available through the package. leaps performs an exhaustive, regression-based search across the proposed variables and selects the model with the highest Adjusted-R2. The algorithm proposes the following 14-factor model with Adjusted-R2 of .9918:

    我已经提出了27个风险因素来解释AHF的回报,但是我不知道是什么构成最佳预测。 可能某些因素无关紧要,从而降低了模型的解释力。 为了选择最佳模型,我使用了可通过软件包使用的基于Adjusted-R2的最佳子集选择算法。 jumps对建议的变量进行详尽的,基于回归的搜索,并选择具有最高Adjusted-R2的模型。 该算法提出以下14因子模型,其中Adj​​usted-R2为.9918:

    • Russell 1000 Value (RLV)

      罗素1000值(RLV)
    • Russell 2000 Index (RUT)

      罗素2000指数(RUT)
    • S&P 500 (GSPC)

      标普500(GSPC)
    • MSCI EAFE (EFA)

      MSCI EAFE(EFA)
    • Barclays Aggregate Bond Index (AGG)

      巴克莱综合债券指数(AGG)
    • Corporate Bond Total Return Index (Corp.TR)

      公司债券总回报指数(Corp.TR)
    • VIX

      VIX
    • Treasury VIX (T.VIX)

      财政部国库券(T.VIX)
    • Convertible Bond Arbitrage Index (CV_ARB)

      可转换债券套利指数(CV_ARB)
    • Equity Market Neutral Index (EQ_Neutral)

      股市中性指数(EQ_Neutral)
    • Multi-Strategy Event Driven Index (MS_EVT)

      多策略事件驱动索引(MS_EVT)
    • Equity Long-Short Index (EQ_LS)

      股票多空指数(EQ_LS)
    • Managed Futures Index (MNGD_FT)

      管理期货指数(MNGD_FT)
    • Multi-Strategy Hedge Fund Index (MS_HF)

      多策略对冲基金指数(MS_HF)

    Now that we have selected our variables, we can estimate the calibrated factor model and see how it does.

    既然我们已经选择了变量,我们就可以估计校准因子模型并查看其效果。

    Image for post

    Based on the results of the regression, we observe that AHF is significantly exposed to traditional sources of risk. Specifically, AHF appears to trade equity and debt and may employ derivatives to either hedge or speculate.

    基于回归结果,我们观察到AHF明显暴露于传统风险源。 具体而言,AHF似乎在买卖股票和债务,并可能使用衍生工具对冲或进行投机。

    Positive exposure to both the S&P 500 (GSPC) and MSCI (EFA) indices suggests that AHF trades global equity and has a long bias. The positive value for AGG further suggests they trade fixed income but may have a slight preference for Treasuries based on the negative coefficient for the corporate bond total return index (corp. tr). The generally significant results for the various hedge fund strategies suggests that AHF employs a complex trading strategy and may use derivatives such as futures (highly significant value for the CS Managed Futures Index, MGND_FT). Futures may be used to either hedge positions or target access to a specific market.

    标普500(GSPC)和MSCI(EFA)指数的正面敞口表明AHF交易全球股票并且具有长期偏见。 AGG的正值进一步表明它们交易固定收益,但基于公司债券总回报指数(corp。tr)的负系数,可能会稍微偏爱国债。 各种对冲基金策略的总体意义重大,这表明AHF采用了复杂的交易策略,并可能使用诸如期货之类的衍生产品(CS管理期货指数MGND_FT的显着价值)。 期货可用于对冲头寸或目标进入特定市场。

    The plot of AHF’s realized returns v. the fitted values from the model demonstrates a high degree of fit and explanatory power.

    AHF的已实现回报与该模型的拟合值的关系图显示出高度的拟合度和解释力。

    Image for post

    C部分:模拟 (Part C: Simulation)

    Parametric and non-parametric Monte Carlo methods are both widely applied in empirical finance, but either presents challenges for estimation.

    参数和非参数蒙特卡罗方法都广泛应用于经验金融中,但都给估计带来挑战。

    Parametric estimation of factor densities requires fitting a large multivariate, fat-tailed probability distribution; which in our specific case would contain 14 variables. Correlations can be notoriously unstable and inaccurate estimation of the variance-covariance matrix would bias the distribution from which we will draw the factor returns. This problem may be overcome by employing copula methods, but this adds to the complexity of the model. On balance, we would prefer to avoid parametric estimation if possible.

    因子密度的参数估计需要拟合较大的多元胖尾概率分布; 在我们的特定情况下,它将包含14个变量。 相关性可能是非常不稳定的,方差-协方差矩阵的不正确估计会偏向于分布我们将得出因子收益的分布。 通过使用copula方法可以解决此问题,但这增加了模型的复杂性。 总而言之,如果可能的话,我们宁愿避免参数估计。

    A potential alternative is non-parametric estimation. To conduct a non-parametric simulation, we could bootstrap the observed, discrete empirical distribution that assigns a probability of 1/T to each of the observed factor returns for t=1…T. This would serve as a proxy to the true density of factor returns and allow us to bypass the messy process of estimating the correlations. However, bootstrap resampling can result in the duplication of some values and the omission of others and while this may be appropriate for inference, it does not provide an obvious advantage in our application.

    潜在的替代方法是非参数估计。 为了进行非参数模拟,我们可以引导观察到的离散经验分布,该分布为t = 1…T的每个观察到的因子收益分配1 / T的概率。 这可以代替要素收益率的真实密度,并允许我们绕过估计相关性的麻烦过程。 但是,自举重采样可能导致某些值重复,而导致其他值遗漏,虽然这可能适合推断,但在我们的应用程序中并未提供明显的优势。

    A more efficient method is simply to take the relatively long history of factor returns as given and add each of the residuals. Simply put, we have 123 months of factor returns (January 2010-March 2020) and 39 residuals (based on the results of the calibration portfolio which spans January 2017-March 2020). If we add each of the 39 residuals to the 123 factor returns, we can produce 123×39 scenarios for the return of AHF (4797 observations in total). This large sample should be capable of providing us with good insight into the tails of AHF’s return distribution and has the advantage of utilizing all of the observed data.

    一种更有效的方法就是简单地考虑给定的相对较长的要素收益历史并添加每个残差。 简而言之,我们有123个月的因子回报(2010年1月至2020年3月)和39个残差(基于跨越2017年1月至2020年3月的校准产品组合的结果)。 如果将39个残差中的每个残差加到123个因子回报中,就可以为123个AHF回报生成123×39个场景(总共4797个观测值)。 这个大样本应该能够为我们提供对AHF收益分布的尾巴的良好洞察,​​并具有利用所有观测数据的优势。

    The simulation proceeds as follows:

    仿真过程如下:

    • Use the calibrated model and factor returns to form predictions for the returns of AHF for January 2010 through March 2020.

      使用校正后的模型和因子收益形成对2010年1月至2020年3月AHF收益的预测。
    • For each of the “i”, i = 1…123, estimated returns, add each of the “j”, j=1…39, residuals.

      对于每个“ i”,i = 1…123,估计收益,将每个“ j”,j = 1…39,残差相加。

    绩效分析 (Performance Analysis)

    Recall that when we introduced this exercise we pretended as though we only had the performance history of AHF from January 2017 through March 2020. Such a short history of performance alone provides only limited insight into the risk/return features of a fund manager over a relatively narrow window of market conditions. To address this shortcoming and provide a more accurate picture of performance we have proposed using Factor Model Monte Carlo (FMMC). The factor model was calibrated using the short, common history of factor and fund returns. The Monte Carlo experiment used factor returns over a longer horizon (January 2010 through March 2020) and the realized factor model residuals to construct 4797 simulated returns for AHF.

    回想一下,当我们引入此练习时,我们假装我们只有从2017年1月到2020年3月才有AHF的业绩历史。仅凭如此短暂的业绩历史,就只能相对有限地了解基金经理的风险/收益特征。市场条件窗口狭窄。 为了解决此缺点并提供更准确的性能描述,我们建议使用因素模型蒙特卡洛(FMMC)。 因子模型是使用较短的常见因子和基金收益历史进行校准的。 蒙特卡洛实验使用更长范围内(2010年1月至2020年3月)的要素收益率和已实现的要素模型残差来构造AHF的4797个模拟收益率。

    To evaluate the performance of our model we will focus on the results for the mean annual return and volatility as well as the venerable Sharpe and Sortino Ratios. Let’s see how we did.

    为了评估模型的性能,我们将重点关注平均年收益率和波动率以及可追溯的夏普和索蒂诺比率的结果。 让我们看看我们是如何做到的。

    1. Average Return

    1.平均回报

    The table below depicts the mean (i.e. average) annual return for the factor model Monte Carlo (FMMC), full history of AHF (January 2010-March 2020) and the truncated/”observed” history (January 2017-March 2020):

    下表描述了因子模型蒙特卡洛(FMMC),AHF的完整历史记录(2010年1月至2020年3月)和截断/“观察到的”历史记录(2017年1月至2020年3月)的平均(即平均)年收益:

    Image for post

    Immediately we can see the improvement that the FMMC model has over the Truncated period. The FMMC model is able to fully capture the return dynamic of AHF while the Truncated return substantially underestimates full history mean.

    马上我们可以看到FMMC模型在截断期间的改进。 FMMC模型能够完全捕获AHF的返回动态,而截断的返回则大大低估了整个历史均值。

    2. Volatility

    2.波动性

    Accurate estimation of the mean alone cannot support the claim that our model is robust. Of equal importance is the volatility. The below table shows the annualized volatility (i.e. standard deviation) for each of the periods under consideration:

    仅凭均值的准确估计就不能支持我们的模型稳健的说法。 波动同样重要。 下表显示了所考虑的每个时期的年度波动率(即标准差):

    Image for post

    Both the FMMC and Truncated estimates slightly undershoot the realized volatility of AHF over the full period. However, both estimated are very close.

    FMMC和“截断”估计都略微低于AHF在整个时期内已实现的波动性。 但是,两者估计都非常接近。

    3. Sharpe Ratio

    3.夏普比率

    With the mean and volatility estimates in hand, we can now calculate the Sharpe Ratio. The Sharpe Ratio is calculated as follows:

    有了平均值和波动率估算值,我们现在可以计算夏普比率。 夏普比率计算如下:

    Image for post

    For most of the test period (January 2010-March 2020) the risk-free rate as proxied by the yield on the 3-month T-Bill was very close to 0%. For simplicity we will adopt 0% as the risk-free rate for our calculations. The below table shows the results:

    在大多数测试期间(2010年1月至2020年3月),由3个月国库券收益率所代表的无风险利率非常接近0%。 为简单起见,我们将采用0%作为无风险利率进行计算。 下表显示了结果:

    Image for post

    The FMMC estimate shows dramatic improvement over the Truncated for estimating the Sharpe Ratio of AHF. This is not necessarily surprising as above we showed the mean return for the Truncated period to be poor while the estimate for the FMMC was quite close. Naturally this will feed into the results for Sharpe, but, again, the results show the utility of the FMMC approach.

    FMMC估计值比截断值(用于估计AHF的Sharpe比率)显着提高。 这不一定是令人惊讶的,因为上面我们显示了截断期的平均收益很低,而FMMC的估算却非常接近。 自然,这将被纳入Sharpe的结果中,但是结果再次显示了FMMC方法的实用性。

    4. Sortino Ratio

    4. Sortino比率

    Finally, we turn to the Sortino Ratio. Sortino is similar to Sharpe, but instead of total volatility, it is focused on what is termed “downside volatility”; or the standard deviation of returns below a stated threshold. Typically, the threshold is set to 0%; the idea being that volatile, positive returns are not considered “bad” because you are still making money, but volatile negative returns suggest an outsized chance of large losses. A higher ratio is considered better. The Sortino Ratio is calculated as follows:

    最后,我们转到Sortino比率。 Sortino与Sharpe类似,但不是总波动率,而是着眼于所谓的“下行波动率”。 或低于规定门槛的收益标准偏差。 通常,阈值设置为0%; 之所以这样的想法是,因为您仍在赚钱,所以波动的正收益不被认为是“坏”的,但是波动的负收益表明出现巨额亏损的可能性很大。 比率越高越好。 Sortino比率计算如下:

    Image for post

    The table depicts the results:

    下表描述了结果:

    Image for post

    The FMMC estimate is very close to the full period and accurately expresses the volatility of the downside returns. We see marked improvement over the Truncated estimate which is lower due to a combination of a lower average return and more downside volatility.

    FMMC的估计非常接近整个周期,并准确表示了下行收益的波动性。 我们看到,由于平均收益较低和下行波动较大,两者的总和比截断后的估计要低得多。

    结论性意见 (Concluding Comments)

    Manager evaluation is one of the oldest and most common problems in investment finance. When the track record of the manager is short it can be difficult to assess the efficacy of the strategy which has ramifications for both fund managers and fund allocators.

    经理评估是投资金融中最古老,最常见的问题之一。 如果经理的业绩记录很短,则可能难以评估该策略的有效性,这对基金经理和基金分配者都有影响。

    In this post, we introduced Factor Model Monte Carlo (FMMC) as a possible solution to the short history problem and used the real world example of Aric’s Hedge Fund (AHF) to demonstrate its efficacy. By using a factor model and the common, short history of fund and factor returns, we estimated the exposure of AHF to different sources of economic and market risk. We were then able to simulate the returns of the AHF to construct a longer history of returns with the goal of gaining improved insight into the fund’s long term performance.

    在本文中,我们介绍了因素模型蒙特卡洛(FMMC)作为短期历史问题的可能解决方案,并使用了Aric对冲基金(AHF)的真实示例来证明其有效性。 通过使用因子模型以及常见的短期资金和因子收益历史,我们估计了AHF暴露于不同的经济和市场风险来源的风险。 然后,我们能够模拟AHF的收益,以构建更长的收益历史,目的是获得对基金长期业绩的更深入了解。

    The results from the FMMC method showed dramatic improvement over using the short history of returns in isolation. Using the full history of returns for AHF as a comparison, we see that the FMMC method accurately models the return, volatility, Sharpe and Sortino Ratios of the fund. By comparison, the truncated history of returns severely underestimated the performance of AHF which would have the consequence of misleading investors.

    FMMC方法的结果表明,与单独使用较短的收益历史相比,有了显着的改进。 使用AHF的全部收益历史作为比较,我们可以看到FMMC方法可以准确地模拟基金的收益,波动率,夏普和Sortino比率。 相比之下,截短的收益历史严重低估了AHF的业绩,这可能会误导投资者。

    Factor Model Monte Carlo has proven to be an effective technique for modeling risk and return for complex strategies and serves as a powerful addition to the fund analyst’s tool kit.

    因子模型蒙特卡洛(Monte Carlo)已被证明是对复杂策略的风险和回报建模的有效技术,并且是基金分析师工具包的有力补充。

    Until next time, thanks for reading!

    直到下一次,感谢您的阅读!

    Aric Lux.

    阿里克斯·勒克斯(Aric Lux)。

    Originally published at http://lightfinance.blog on August 28, 2020.

    最初于 2020年8月28日 发布在 http://lightfinance.blog 上。

    翻译自: https://towardsdatascience.com/better-portfolio-performance-with-factor-model-monte-carlo-in-r-3910d0a6ceb6

    风险评估模型蒙特卡洛模型

    展开全文
  • 发现物质功率谱的增强是这些模型中的通用特征。 特别地,我们表明,在前一种类型(例如Starobinsky模型)中,光谱的放大倍数比后一种(例如指数模型)大得多。 在可行的f(R)模型中,允许的中微子总质量(m)的...
  • datamodelr:R中的数据模型
  • 介绍 聚类模型是一个概念,用于表示我们试图识别聚类类型。四种最常见聚类方法模型是层次...基于模型的聚类是迭代方法,通过优化聚类数据集分布,将一组数据集拟合到聚类。高斯分布只不过是正态分...

    原文链接:http://tecdat.cn/?p=6105

    原文出处:拓端数据部落公众号

     

    介绍

     

    四种最常见的聚类方法模型是层次聚类,k均值聚类,基于模型的聚类和基于密度的聚类 

    可以基于两个主要目标评估良好的聚类算法:

    • 高组内相似性
    • 低组间相似性

     

    基于模型的聚类是迭代方法,通过优化聚类中数据集的分布,将一组数据集拟合到聚类中。高斯分布只不过是正态分布。此方法分三步进行:

    1. 首先随机选择高斯参数并将其拟合到数据集。
    2. 迭代地优化分布参数以拟合尽可能多的点。
    3. 一旦收敛到局部最小值,您就可以将数据点分配到更接近该群集的分布。

    高斯混合模型

    基于概率模型的聚类技术已被广泛使用,并且已经在许多应用中显示出有希望的结果,从图像分割,手写识别,文档聚类,主题建模到信息检索。基于模型的聚类方法尝试使用概率方法优化观察数据与某些数学模型之间的拟合。

    生成模型通常使用EM方法求解,EM方法是用于估计有限混合概率密度的参数的最广泛使用的方法。基于模型的聚类框架提供了处理此方法中的几个问题的主要方法,例如密度(或聚类)的数量,参数的初始值(EM算法需要初始参数值才能开始),以及分量密度的分布(例如,高斯分布)。EM以随机或启发式初始化开始,然后迭代地使用两个步骤来解决计算中的循环:

    • E-Step。使用当前模型参数确定将数据点分配给集群的预期概率。
    • M-Step。通过使用分配概率作为权重来确定每个集群的最佳模型参数。

     

    R中的建模

     

    mb = Mclust(iris[,-5])
    
    #或指定集群数
    mb3 = Mclust(iris[,-5], 3)
    
    # 或指定集群数
    mb$modelName
    
    # 或指定集群数
    mb$G
    
    # 在给定聚类中观察的概率
    head(mb$z)
    
    # 获取概率,均值,方差
    summary(mb, parameters = TRUE)
    

     

    table(iris$Species, mb$classification)
    # 比较
    table(iris$Species, mb3$classification)
    

    比较每个群集中的数据量

    在将数据拟合到模型中之后,我们基于聚类结果绘制模型。

    让我们绘制估计的密度。 

    plot(mb, "density")
    

     

    您还可以使用该summary()函数来获取最可能的模型和最可能数量的集群。对于此示例,最可能的簇数为5,BIC值等于-556.1142。

    比较聚类方法

    在使用不同的聚类方法将数据拟合到聚类中之后,您可能希望测量聚类的准确性。在大多数情况下,您可以使用集群内集群度量标准作为度量。集群间距离越高越好,集群内距离越低、越好。

    接下来,检索聚类方法的集群验证统计信息:

    通常,我们使用within.cluster.ssavg.silwidth验证聚类方法。该within.cluster.ss测量表示所述簇内总和的平方,和avg.silwidth表示平均轮廓宽度

    • within.cluster.ss测量显示了相关对象在群集中的紧密程度; 值越小,集群中的对象越紧密。
    • avg.silwidth是一种度量,它考虑了群集中相关对象的紧密程度以及群集之间的区别程度。轮廓值通常为0到1; 接近1的值表明数据更好地聚类。

     

    k-means和GMM之间的关系

    K均值可以表示为高斯混合模型的特例。通常,高斯混合有更好的表现力,因为数据的群集分配还取决于该群集的形状,而不仅仅取决于其接近度。

    与k-means一样,用EM训练高斯混合模型可能对冷启动条件非常敏感。如果我们将GMM与k-means进行比较和对比,我们会发现前者的初始条件比后者更多。 

    结果

     每个聚类被建模为多元高斯分布,并通过给出以下内容来指定模型:

    1. 集群数量。
    2. 每个群集中所有数据点的分数。
    3. 每个聚类的均值和它的协方差矩阵。

      

    参考文献

    1.R语言k-Shape算法股票价格时间序列聚类

    2.R语言中不同类型的聚类方法比较

    3.R语言对用电负荷时间序列数据进行K-medoids聚类建模和GAM回归

    4.r语言鸢尾花iris数据集的层次聚类

    5.Python Monte Carlo K-Means聚类实战

    6.用R进行网站评论文本挖掘聚类

    7.用于NLP的Python:使用Keras的多标签文本LSTM神经网络

    8.R语言对MNIST数据集分析 探索手写数字分类数据

    9.R语言基于Keras的小数据集深度学习图像分类

    展开全文
  • MMFAIR-R中的农业混合模型 一方面可以被视为更大的兄弟 ,因为它处理更复杂的分析。 但是,另一方面,它也解决了R中(广义)线性混合模型分析的局限性: 迄今为止,R还没有用于此类分析的单独入门套件。 相反, ...
  • mlrMBO:R中的贝叶斯优化和基于模型的优化的工具箱
  • R中模型性能提升

    千次阅读 2016-04-29 18:00:34
    最近正在学习这一部分,备忘录。1、参数调整 对算法合适选项进行调整过程...(1)训练哪种模型,(2)模型中哪些参数可调,可调节空间多大,(3)选择评价标准 以C5.0示例:library(caret)control <- trainCont

    最近正在学习这一部分,备忘录。

    1、参数调整
    对算法合适的选项进行调整的过程——参数调整
    caret包中提供了多种工具进行自动参数调整,train()函数作为接口,可以选择评估方法和度量性指标,自动寻优过程。
    主要考虑的问题:
    (1)训练哪种模型,(2)模型中哪些参数可调,可调节空间多大,(3)选择评价标准
    以C5.0示例:

    library(caret)
    
    control <- trainControl(method = "cv",number =5,
                            selectionFunction = "oneSE")
    set.seed(7)
    fit <-train(label~.,data = traindata1[-1],method ="C5.0",
                metric = "Accuracy",trControl = control)
    print(fit)

    定制调整的过程,trainControl()函数实现配置选项,主要参数为method,selectionFunction。?trainControl可以得到所有的参数列表。method主要用来设置重抽样的方法,例如保持抽样或k折验证。selectionFunction可设定一个函数在各个候选者中选择最优的模型。
    可创建可调参数的数据框,每一行代表可调参数的组合

    grid <- expand.grid(.model = “tree”,
    .trials = c(1,5,10,15,20,25,30,35),
    .winnow = “FALSE”
    )
    set.seed(300)
    m <- train(label~.,data = traindata1[-1],method = “C5.0”,
    metric = “Kappa”,
    trControl = ctrl,
    tuneGrid = grid)

    2、使用元学习提高模型性能
    bagging: 自主抽样产生数据集,这些数据集使用单一机器学习算法产生多个模型,最后采用投票或平均的方法来组合预测值
    boosting:在不同的重抽样数据中训练模型的集成学习。区别:重抽样数据集是专门构建的,根据上一次分类器,自动调整本次分类器构造时,各样本被抽中的概率。对经常分错的样本给与较大的权重,预测正确的样本在下一个分类器中的可能性较小。通过投票决定最后的预测值,每个分类器的票数会按照建模数据集上的准确度进行加权。

    以下主要介绍bagging:
    ipred包提供了bagging决策树的实现。如下:

    library(ipred)

    m.bag <- bagging(label~.,data = traindata1[,-1], nbagg = 25)
    p.bag <- predict(m.bag,testdata1[,-c(1,16)],type = “class”)

    其中nbagg参数控制用来投票的决策树的数目。
    caret包中train函数采用10折交叉验证方法建立bagging树,函数为treebag

    library(caret)

    control <- trainControl(method = “cv”,number =10)

    set.seed(7)
    fit <-train(label~.,data = traindata1[-1],method =”treebag”,
    metric = “Accuracy”,trControl = control)

    caret包中含添加更加通用的bag()函数,支持向量机svmBag,朴素贝叶斯nbBag,决策树ctreeBag,神经网络nnetBag。以下选取支持向量机为例svmBag,创建一个bagging控制对象:

    ctrl <- trainControl(method = “cv”,number = 10)
    bagctrl <- bagControl(fit = svmBagfit,predict=svmBagpred,
    aggregate = svmBagaggregate)set.seed(300)svmbag<train(label .,data=traindata1[,1],bag,trControl=ctrl,bagControl=bagctrl)使baggingsvmBagfit
    function (x, y, …)
    {
    loadNamespace(“kernlab”)
    out <- kernlab::ksvm(as.matrix(x), y, prob.model = is.factor(y),
    …)
    out
    }

    以及svmBagpred,svmBagaggregate的函数实现方式,可定制bagging

    展开全文
  • R语言中的模型修正函数update

    千次阅读 2016-03-20 22:17:01
    R语言中的模型修正函数update的使用方法
  • rstudio 创建模型Creating a model is an essential part of forecasting and data analysis. I’ve put together a quick guide on my process for ...创建模型是预测和数据分析重要组成部分。 我已经为我数...
  • 在实践, 因子负载较低(或测量质量较差)的模型的拟合指数要好于因子负载较高的模型。例如,如果两个模型具有相同错误指定级别,并且因子负载为.9的模型的RMSEA可能高于.2,而因子负载为.4的模型的RMSEA可能...
  • 预测模型 使用时间序列模型R中的英国GDP进行预测 使用的模型:ARIMA,auto.arima,Naive,ETS 对于模型性能评估,考虑了Diebold / Mariano测试和RMSE。 上传了项目摘要doc文件,以供详细参考。
  • 使用R语言中的GWmodel进行GWR模型的运算 一、安装R语言 1、R安装包下载地址:https://cran.r-project.org/ 2、RStudioa安装包下载地址:https://rstudio.com/products/rstudio/ 二、使用步骤 1.安装GWmodel包 打开...
  • R语言混合线性模型的实现以及参数解析

    万次阅读 热门讨论 2019-04-10 20:09:39
    前言 为什么要用混合线性模型:比如测量了不同收入水平人群收入和幸福感,但每个群体内收入水平是不同,幸福感也不同,两者之间关系也是不同, 如果直接用一般线性模型,会造成错误...R中混合线性模型...
  • R语言中的模型公式与图表

    千次阅读 2015-01-19 15:15:57
    #利用base包cars数据集来讲解模型和公式这一章 cars.lm(formula=dist~speed,data=cars) cars.lm #其中lm是估算线性模型参数函数,formula代表公式,dist~speed就是公式,cars是数据集 #使用summary函数 summary...
  • 这些练习对我来说尤其具有启发性,因为它们说明包含随机效应(又称变化效应)不仅可以改变相对模型排名,而且还强调,添加随机效应可以极大地改变我们对固定效应的估计(即,通常情况下,关心我们的模型中的大多数...
  • 建立流失玩家预测模型之初,对于P值和F值不是非常理解,后来随着模型的建立,清楚了P值和F值意义,结合实际业务,对这两个值进行权衡。 P值,英文为precision,准确率。 R值,英文为recall,召回率。 实例:...
  • R语言进行期权定价Heston模型

    千次阅读 2020-03-02 17:47:02
    在本文,我将向您展示如何模拟股票价格Heston随机波动率模型。 Heston模型是是一种期权估值方法,它考虑到同一资产在给定时间交易不同期权波动性变化。它试图通过使用随机过程来模拟波动率和利率来重新创建...
  • R语言中的划分聚类模型

    千次阅读 2019-09-09 17:30:35
    划分聚类是用于基于数据集相似性将数据集分类为多个组聚类方法。 分区聚类,包括: K均值聚类(MacQueen 1967),其中每个聚类由属于聚类数据点中心或平均值表示。K-means方法对异常数据点和异常值敏感。 ...
  • rugarch包与R语言中的garch族模型

    千次阅读 2018-03-24 17:54:00
    rugarch包是R中用来拟合和检验garch模型的一个包。该包最早在http://rgarch.r-forge.r-project.org上发布,现已发布到CRAN上。简单而言,该包主要包括四个功能: 拟合garch族模型 garch族模型诊断 garch族模型...
  • R语言系统教程(十二):R语言与线性模型有关函数12.1 基本函数12.2 提取模型信息通用函数 12.1 基本函数 适应于多元线性模型的基本函数是lm(),其调用形式是 lm(formula, data = data.frame) 其中formula为...
  • 概率图模型 基于R语言 这本书中的第一个R语言程序 1 prior <- c(working =0.99,broken =0.01) 2 likelihood <- rbind(working = c(good=0.99,bad=0.01),broken =c(good=0.6,bad=0.4)) 3 data <- c...
  • R:怎么在混合模型中分析随机效应显著性 对于农业试验研究,试验设计经常是随机区组(Randomized Block Design)或裂区设计(Split Plot Design),在这些试验基础上产生数据进行方差分析需要用到混合模型(Mixed...
  • R中的几种统计分布及常用模型

    千次阅读 2018-05-04 10:26:49
    简单介绍一下R中的几种统计分布及常用模型 转载自:https://www.cnblogs.com/nxld/p/6060360.html 统计学上分布有很多,在R中基本都有描述。因能力有限,我们就挑选几个常用的、比较重要的简单介绍一下每种分布的...
  • R语言中的隐马尔可夫HMM模型实例

    千次阅读 2020-04-10 22:39:53
    HMM用于建模数据序列,无论是从连续概率分布还是从离散概率分布得出。它们与状态空间和高斯混合模型相关,因为它们旨在估计引起观测状态。状态是未知或“隐藏”,并且HMM试图估计状态,类似于无监督聚类过程。
  • 在本教程,您将学习如何在R中创建神经网络模型。 神经网络(或人工神经网络)具有通过样本进行学习能力。人工神经网络是一种受生物神经元系统启发信息处理模型。它由大量高度互连处理元件(称为神经元)...
  • 本文 将针对R进行几次建模练习结果,以魁北克数据为依据,分为13年训练和1年测试。prophet与基本线性模型(lm),一般加性模型(gam)和随机森林(randomForest)进行了比较。 首先,设置一些选项,加载...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 12,667
精华内容 5,066
关键字:

模型中的r