0 library

library(ggplot2)
library(gridExtra)
library(stats)

1 MA(1)模型的ACVF,ACF,PACF

1.1生成高斯白噪声

set.seed("1234")
n=250
a=rnorm(n)

1.2生成含250个值的序列,分布如图

u=0
Theta=0.5
z=array(NA,n)
z[1]=u+a[1]
for(t in 2:n) 
  z[t]=u+a[t]-Theta*a[t-1]
plot(z,type = "l")

1.3手工计算ACVF,ACF估计值

Ez=mean(z)
Gamma=array(NA,n-1)
Rou=array(NA,n-1)
sum=0
for(t in 1:n) 
  sum=sum+(z[t]-Ez)*(z[t]-Ez)
Gamma0=1/n*sum

for(k in 1:(n-1)){
  sum=0
  for(t in 1:(n-k)) 
    sum=sum+(z[t]-Ez)*(z[t+k]-Ez)
  Gamma[k]=1/n*sum
  Rou[k]=Gamma[k]/Gamma0
}                                      

No=0:(n-1)
data=data.frame(no = No, gamma=c(Gamma0,Gamma),rou = c(1,Rou))

p <-ggplot(data) 
p1=p+geom_bar(width=0.1,mapping=aes(x=no, y=rou),stat = "identity")+
  xlab("NO")+ ylab("ACF")+ xlim(-1,24)+geom_hline(yintercept = 0)
p2=p+geom_bar(width=0.1,mapping = aes(x=no, y=gamma),stat ="identity")+
  xlab("NO")+ ylab("ACVF")+ xlim(-1,24)+geom_hline(yintercept = 0)
grid.arrange(p1, p2, nrow=1, ncol=2)

1.4 acf函数生成

acf(z,type = "correlation")

acf(z,type = "covariance")

acf(z,type = "partial")

2.1 AR(1)模型的标准谱密度函数

f <- function(w,b) (1/(2*pi))*((1-b^2)/(1+b^2-2*b*cos(w)))

# 创建数据点
w <- seq(0, pi, by=0.01)
y1 <- f(w,0.7)
y2 <- f(w,0.5)
y3 <- f(w,0.3)
df <- data.frame(w, y1, y2, y3)

# 用ggplot2来画图
g <- ggplot(df)
g1 <- g + geom_line(col='red',mapping=aes(w,y1),linetype = 1)+
  geom_line(col='blue',mapping=aes(w,y2),linetype = 2)+
  geom_line(col='black',mapping=aes(w,y3),linetype = 3)+
  geom_hline(yintercept = 0)+geom_vline(xintercept = 0)+
  ylab("f(w)")+
  ggtitle("f(w)")
print(g1)

2.2谱密度估计

三种常用方法:

1.(最小二乘)平滑周期图

    对有限样本,异方差性影响平滑的有效性

2.(最小二乘)平滑对数周期图

3.局部似然估计(推荐这个,因为前两个不是有效的估计)
## Examples from Venables & Ripley
## spec.pgram 

##Allowed methods are "pgram" (the default) and "ar".
##这是啥方法

par(mfrow = c(2,2))

spectrum(lh)
spectrum(lh, spans = 3)
spectrum(lh, spans = c(3,3))
spectrum(lh, spans = c(3,5))

spectrum(ldeaths)
spectrum(ldeaths, spans = c(3,3))
spectrum(ldeaths, spans = c(3,5))
spectrum(ldeaths, spans = c(5,7), log = "dB", ci = 0.8)

# for multivariate examples see the help for spec.pgram

## spec.ar
spectrum(lh, method = "ar")
spectrum(ldeaths, method = "ar")

关于R新增的了解

1.```{r,warning=FALSE,dev=‘svg’} 输出矢量图

2.ggplot2 的使用

3.使用ggplot2时,常用的par不能用了,可以使用gridExtra包里的grid.arrage分页

4.acf,spectrum函数的使用

5.ggplot2的坐标轴画法

6.set.seed(“1234”) 设置伪随机数起点,使每次模拟结果相同

另外:结构断点与谱分析http://www.xzbu.com/3/view-4489847.htm