用R语言软件估计光谱密度
任何时间序列都可以表示为在基波(谐波)频率上振荡的余弦和正弦波之和= j / n,其中j = 1,2,...,n / 2。周期图给出了关于各种频率的相对强度的信息,用于解释时间序列的变化。
周期图是称为谱密度的群体函数的样本估计,其是群体平稳时间序列的频域表征。谱密度是与自协方差时域表示直接相关的时间序列的频域表示。本质上,谱密度和自协方差函数包含相同的信息,但以不同的方式表达。
[ 回顾注释:自协方差是自相关的分子。自相关是自协方差除以方差。]
假设γ(h)是静止过程的自协方差函数,f(ω)是同一过程的谱密度。在前一句的符号中,h =时滞和ω=频率。
在高级微积分的语言中,自协方差和谱密度是傅立叶变换对。我们不会担心这种情况的微积分。我们将重点关注谱密度的估计 - 一系列的频域特征。这里仅给出傅里叶变换方程以确定在时域表示和序列的频域表示之间存在直接链接。
在数学上,频谱密度是针对负频率和正频率定义的。然而,由于函数的对称性及其对于-1/2到+1/2范围之外的频率的重复模式,我们只需关注0和+1/2之间的频率。
“总”积分谱密度等于系列的方差。因此,特定频率区间内的谱密度可以被视为由那些频率解释的方差的量。
估计光谱密度的方法
原始周期图是人口谱密度的粗略样本估计。估计是“粗略的”,部分是因为我们只使用离散的基波谐波频率用于周期图,而频谱密度是在连续的频率上定义的。
对谱密度的周期图估计的一种可能的改进是使用居中移动平均值来平滑它。可以使用逐渐减小的方法创建另外的“平滑”,该方法对系列的末端(时间)加权小于数据的中心。我们不会在本课中介绍逐渐减少的内容。有兴趣的人可以参阅本书第4.4节和各种互联网资源。
平滑周期图的另一种方法是基于以下事实的参数估计方法:任何静止时间序列可以通过某种顺序的AR模型来近似(尽管它可能是高阶)。在该方法中,找到合适的AR模型,然后将谱密度估计为该估计的AR模型的谱密度。
平滑方法(光谱密度的非参数估计)
平滑周期图的常用方法具有如此精美的名称,听起来很难。实际上,它只是一个集中的移动平均程序,只有一些可能的修改。对于时间序列,具有参数m的Daniell核是居中移动平均值,其通过平均时间t-m和t + m(包括)之间的所有值在时间t创建平滑值。例如,m = 2 的Daniell内核的平滑公式为
^ x t= x t - 2 + x t - 1 + x t + x t + 1 + x t + 2 5x^t=xt−2+xt−1+xt+xt+1+xt+25
在R中,可以使用命令kernel(“daniell”,2)生成m = 2 的Daniell内核的加权系数。 结果是
coef [-2] = 0.2
coef [-1] = 0.2
coef [0] = 0.2
coef [1] = 0.2
coef [2] = 0.2
coef []的下标是指与时间t的平均值中心的时差。因此,在这种情况下的平滑公式是
^ x t=0.2xt-2+0.2xt-1+0.2xt+0.2xt+1+0.2xt+2,x^t=0.2xt−2+0.2xt−1+0.2xt+0.2xt+1+0.2xt+2,
这与上面给出的公式相同。
修改后的Daniell内核使得平均值中的两个端点接收内部点的重量的一半。对于m = 2 的修改后的Daniell内核,平滑是
^ x t= x t - 2 +2 x t - 1 +2 x t +2 x t + 1 + x t + 2 8= 0.125 x t - 2 + 0.25 x t - 1 + 0.25 x t + 0.25 x t + 1 + 0.125 x t + 2x^t=xt−2+2xt−1+2xt+2xt+1+xt+28=0.125xt−2+0.25xt−1+0.25xt+0.25xt+1+0.125xt+2
在R中,命令内核(“modified.daniell”,2)将列出刚刚使用的加权系数。
可以对Daniell内核或修改后的Daniell内核进行卷积(重复),以便将平滑再次应用于平滑值。通过在更宽的时间间隔内取平均值,可以产生更广泛的平滑效果。例如,为了重复丹尼尔内核米上起因于一个丹尼尔内核与平滑值= 2 米 = 2时,公式将是
^ ^ x t= ^ x t - 2 + ^ x t - 1 + ^ x t + ^ x t + 1 + ^ x t + 2 5x^^t=x^t−2+x^t−1+x^t+x^t+1+x^t+25
这是在任一方向上的两个时间段t内的平滑值的平均值。
在R中,命令内核(“daniell”,c(2,2))将提供系数,这些系数将作为权重应用于对两个平滑中m = 2 的回旋Daniell内核的原始数据值求平均值。结果是
coef [-4] = 0.04
coef [-3] = 0.08
coef [-2] = 0.12
coef [-1] = 0.16
coef [0] = 0.20
coef [ 1] = 0.16
coef [2] = 0.12
coef [3] = 0.08
coef [4] = 0.04
这会生成平滑公式
^ x t=0.04xt-4+0.08xt-3+0.12xt-2+0.16xt-1+0.20xt+0.16xt+1+0.12xt+2+0.08xt+3+0.04xt+4。x^t=0.04xt−4+0.08xt−3+0.12xt−2+0.16xt−1+0.20xt+0.16xt+1+0.12xt+2+0.08xt+3+0.04xt+4.
其中端点具有较小重量的修改方法的卷积也是可能的。命令内核(“modified.daniell”,c(2,2)) 给出了这些系数:
coef [-4] = 0.01563
coef [-3] = 0.06250
coef [-2] = 0.12500
coef [-1] = 0.18750
coef [0] = 0.21875
coef [1] = 0.18750
coef [2] = 0.12500
coef [3] = 0.06250
coef [4] = 0.01563
因此,中心值的加权比未修改的Daniell内核略重。
当我们平滑周期图时,我们在频率间隔而不是时间间隔上进行平滑。请记住,周期图是在基本频率ω确定Ĵ = Ĵ/ N为Ĵ = 1,2,...,ñ / 2。让我(ω Ĵ )表示在频率ω的周期图值Ĵ = Ĵ/ N。当我们使用带参数m的Daniell内核来平滑周期图时,平滑值\(\ hat {I}(\ omega_j)\)是范围(jm)/ n到(j)中频率的周期图值的加权平均值+ m)/ n。^ 我(ωĴ)I^(ωj)
带宽
带宽应该足以平滑我们的估计,但是如果我们使用太大的带宽,我们将过多地平滑周期图并且错过看到重要的峰值。在实践中,通常需要一些实验来找到提供合适平滑的带宽。
带宽主要由平滑中平均值的数量控制。换句话说,Daniell内核的m参数以及内核是否被卷积(重复)会影响带宽。
注意: 带有图表的带宽R报告与使用上述公式计算的值不匹配。请参阅第12页的脚注。190您的文字作为解释。
R代码
使用Daniell内核对周期图进行平均/平滑可以使用两个命令的序列在R中完成。第一个定义了Daniell内核,第二个定义了平滑周期图。
例如,假设观察到的序列被命名为x,我们希望使用具有m = 4的Daniell内核来平滑周期图。命令是
mvspec(x,k,log =“no”)
第一个命令创建平滑所需的加权系数,并将它们存储在名为k的向量中。(将它称为k是任意的。它可以被称为任何东西。)第二个命令要求基于系列x的周期图的谱密度估计,使用存储在k中的加权系数,该图将是普通的比例,不是对数刻度。
如果需要卷积,则可以将内核命令修改为类似k = kernel(“daniell”,c(4,4))的内容。
有两种方法可以实现修改后的Daniell内核。您可以更改内核命令以引用“modified.daniell”而不是“daniell”,也可以跳过使用内核命令并在mvspec命令中使用spans参数。
spans参数给出了所需修改的Daniell内核的长度(= 2 m +1)。例如,m = 4 的修改后的Daniell内核长度L = 2 m +1 = 9,因此我们可以使用该命令
mvspec(x,spans = 9,log =“no”)可以使用每次通过m = 4 的修改的Daniell内核的两次通过mvspec(x,spans = c(9,9),log =“no”)
示例:此示例将使用文本中多个位置使用的鱼类招募系列,包括第4章中的几个位置。该系列包含n = 453个月度值,用于衡量南半球位置的鱼群数量。数据位于文件recruit.dat中。
可以使用命令创建原始周期图(或者可以使用第6课中给出的方法创建原始周期图)。
mvspec(x,log =“no”)
请注意,在刚刚给出的命令中,我们省略了为平滑提供权重的参数。
原始周期图如下:
下一个图是使用具有m = 4 的Daniell核的平滑周期图。注意,平滑的一个效果是未平滑版本中的主峰现在是第二高峰。之所以发生这种情况是因为峰值在未平滑的版本中如此清晰地定义,当我们用几个周围值对其进行平均时,高度会降低。
下一个图是一个平滑的周期图,使用Daniell内核的两次传递,每次传递m = 4。请注意它是如何比以前更平滑。
要了解两个主要峰的位置,请为mvspec输出指定一个名称,然后列出它。例如,
specvalues = mvspec(x,k,log =“no”)
specvalues $ spec
您可以筛选输出以查找峰值出现的频率。频率和频谱密度估计单独列出,但顺序相同。确定最大光谱密度,然后找到相应的频率。
这里,第一个峰值的频率≈0208。与此周期相关的周期(月数)= 1 / .0208 = 48个月。第二个峰值出现在频率≈0.083333。相关时期= 1 / .08333 = 12个月。第一个峰值与厄尔尼诺天气效应有关。第二个是通常12个月的季节性影响。
这两个命令将垂直虚线放在峰值密度的近似位置处的(估计的)谱密度图上。
abline(v = 1/44,lty =“dotted”)
abline(v = 1/12,lty =“dotted”)
这是结果图:
我们已经足够平滑,但出于演示目的,下一个情节是结果
mvspec(x,spans = c(13,13),log =“no”)
这使用两次修改的Daniell内核,每次长度L = 13(所以m = 6)。情节有点平滑,但不是很多。顺便说一句,峰值与上面的图中完全相同。
绝对有可能过度平滑。假设我们使用总长度= 73(m = 36)的修改过的Daniell内核。命令是
mvspec(x,spans = 73,log =“no”)
结果如下。高峰消失了!
光谱密度的参数估计
谱密度估计的平滑方法称为非参数方法,因为它不使用任何参数模型用于基础时间序列过程。另一种方法是参数方法,该方法需要找到该系列的最佳拟合AR模型,然后绘制该模型的谱密度。该方法由一个定理支持,该定理表明任何时间序列过程的谱密度可以通过AR模型的谱密度(某种顺序,可能是高顺序)来近似。
在R中,使用命令/功能spec.ar可以轻松完成谱密度的参数估计。像spec.ar(x,log =“no”)这样的命令将导致R完成所有工作。同样,为了识别峰值,我们可以通过执行specvalues = spec.ar(x,log =“no”)之类的操作为spec.ar结果指定名称。
对于鱼类补充实例,结果如下。注意,绘制的密度是AR(13)模型的密度。我们当然可以为这些数据找到更简约的ARIMA模型。我们只是使用该模型的谱密度来近似观察到的序列的谱密度。
估计的谱密度的出现与以前相同。估计的厄尔尼诺峰位于略微不同的位置 - 频率约为0.024,周期约为1 / 0.024 =约42个月。
德趋势
在进行光谱分析之前,应该对一系列进行去趋势分析。趋势将导致在低频率下的这种主要谱密度,从而不会看到其他峰值。默认情况下,R命令mvspec使用线性趋势模型执行去趋势。也就是说,使用来自回归的残差来估计谱密度,其中y变量=观测数据和x变量= t。如果存在不同类型的趋势,例如二次方,则可以使用多项式回归来在探索估计的谱密度之前使数据去趋势化。需要注意的是R命令spec.ar,但默认情况下不进行解趋势。
平滑器在原始数据中的应用
请注意,此处描述的平滑器也可以应用于原始数据。Daniell内核及其修改只是移动平均(或加权移动平均)平滑器。
还有问题吗?联系我们!
大数据部落 -中国专业的第三方数据服务提供商,提供定制化的一站式数据挖掘和统计分析咨询服务
统计分析和数据挖掘咨询服务:y0.cn/teradat(咨询服务请联系官网客服)
QQ:3025393450
【服务场景】
科研项目; 公司项目外包;线上线下一对一培训;数据采集;学术研究;报告撰写;市场调查。
【大数据部落】提供定制化的一站式数据挖掘和统计分析咨询服务
转载于:https://www.cnblogs.com/tecdat/p/10524347.html
用R语言软件估计光谱密度相关推荐
- r统计建模与r软件期末考试题_“统计学诺贝尔奖”授予 R 语言软件工程师 Hadley Wickham | 科研圈日报...
"科研圈日报"主要关注科研圈与研究者个体.科研圈与更广阔的社会环境之间的重要互动.点击 这里 可以查看往期内容. · 学术荣誉 "统计学诺贝尔奖"授予 R 语言 ...
- R 语言怎么保存工作目录到当前路径_【R语言基础】01.R语言软件环境搭建及常用操作...
一.R语言简介 R语言是专业的统计分析软件,来自著名数据科学网站(http://www.kdnuggets.com/)发起的一个2019年统计分析和数据挖掘软件使用情况的调查结果: 表明R语言是该领域 ...
- r语言软件GDINA_finTech MSc代做、代写Python程序语言、代写MSc program、代做Python设计帮做C/C++编程|代写R语言...
finTech MSc代做.代写Python程序语言.代写MSc program.代做Python设计帮做C/C++编程|代写R语言Strathclyde Business School, finTe ...
- r语言软件GDINA_认知诊断分析系统(flexCDMs)设计及其实现
[Abstract] Cognitive diagnostic assessment arose in 80 s last century as a new testing mode which co ...
- 云服务器linux(centos)系统下Rstudio的下载及连接R语言软件
前提为在云服务器端打开Rstudio8787的端口 1 进入Rstudio官网找到对应系统的Rstudio的rpm文件用wget下载RStudio Server - Posit好吧Rstudio官网改 ...
- geeglm函数R语言广义估计方程系数的置信区间计算
转载 链接: https://www.codesd.com/item/confidence-interval-of-coefficients-using-the-generalized-estimat ...
- R语言随机搜索变量选择SSVS估计贝叶斯向量自回归(BVAR)模型
介绍 最近我们被客户要求撰写关于向量自回归的研究报告,包括一些图形和统计输出.向量自回归(VAR)模型的一般缺点是,估计系数的数量与滞后的数量成比例地增加.因此,随着滞后次数的增加,每个参数可用的信息 ...
- Day01零基础自学R语言(最详细教程)——R软件安装
R语言简介 R语言是当今排名进入前十五的程序设计语言,也是大数据处理的常用工具之一. R语言是由新西兰奥克兰大学的Ross Ihaka和Robert Gentleman所开发的,因为两人名字开头都是R ...
- 【视频】什么是Bootstrap自抽样及应用R语言线性回归预测置信区间实例|数据分享
最近我们被客户要求撰写关于Bootstrap的研究报告,包括一些图形和统计输出. 自抽样统计是什么以及为什么使用它? 本文将自抽样方法与传统方法进行比较,并了解它为何有用.并在R语言软件中通过对汽车速 ...
- R语言实现向量自回归VAR模型
澳大利亚在2008 - 2009年全球金融危机期间发生了这种情况.政府发布了一揽子刺激计划,其中包括2008年12月的现金支付,恰逢圣诞节支出.因此,零售商报告销售强劲,经济受到刺激,收入增加了. 最 ...
最新文章
- 用keil烧写现成的hex文件
- Python eval()函数的使用
- XP的DNS服务器(BIND)配置
- 2019数据安装勾选_万能的XY数据标签插件,柱形图也可以呈现变化率
- 让表单文本框只读不可编辑的方法
- Java并发编程:如何创建线程?
- 我国快递年业务量首次突破千亿件大关
- volatile指令重排_面试:为了进阿里,重新翻阅了Volatile与Synchro
- 时间序列分析导论书摘:自相关图意义分析
- 【云周刊】第200期:云栖专辑 | 阿里开发者们的第6个感悟:享受折磨
- C# 面向对象初级 (参考传智播客视频)
- 陈国君Java程序设计基础笔记和习题
- 四川山海蓝图抖音上热门的技巧
- C# - 音乐小闹钟_BetaV1.0
- 自然语言处理(NLP)的一般处理流程!
- OpenGL学习笔记:光照贴图
- 用ADC0809实现八通道采集
- 购买汽车都有哪些费用,以及计算公式
- 中国SAP顾问在美国的跳槽经历
- IKEv2的认证数据生成过程
热门文章
- 安卓模拟经营类游戏_十大最诱人手机模拟经营类游戏专题
- h2o api java_H2O框架简介
- zigbee协议栈初使用(四)无线串口透传
- 美国佐治亚理工计算机专业,世界大学排名之:美国佐治亚理工学院
- 操作系统虚拟存储管理实验
- Secondary NameNode工作原理
- 学而思pythonlevel3_【学而思网校语言学习】学而思网校【2019-寒】AE英语直播班 Level 3上【报价 价格 评测 怎么样】 -什么值得买...
- VB如何只读取字符串中的数字部分??
- Javascript中大于和小于
- 李想的理想,不太「理想」