3-代码

#布朗运动
library(somebm)
library(readxl)
library(TSA)
library(forecast)
library(tseries)
library(fUnitRoots)

xlsx_1<-read_excel("~/Desktop/graduation design/WIND/export for USA.xlsx", sheet = "Sheet2")
mon <- xlsx_1[,2]
mon.timeseries <- ts(mon,start = c(1995,1),end = c(2019,12),frequency = 12)
rt = diff(log(mon.timeseries))    #299个数据
#计算残差
aipx = rep(0,299)

for (i in 13:299){
  aipx[i] = rt[i]+0.0366*rt[i-1]+0.5681*aipx[i-1]+0.4999*aipx[i-12] - cos(i/24);
}


#构造CUSUM检测函数
D <- function(k){ ## function关键字
  res = 0
  for(i in 1:k){
    res = res + aipx[60+i]
  }
  for(i in 1:299){
    res = res - k/299 * aipx[i]
  }
  return(res)       ## 返回值
}

#gam = 0的阈值函数
g <- function(k){
  return((1+k/m)*(k/m + k)^0)
}

T = 3               #检测集与训练集的比例 
m = 60              #训练集个数
c = 0.05            #显著性水平
k = m * T           #检测集个数
tao = rep(0,k)      #停时向量
rho = 0.006972      #方差
eta = aipx^2-1      
eta2 = mean(eta^2)
for(i in 1:k){
  j = 1
  if(abs(D(j))/(sqrt(m) * rho * g(j) >= c)){
    tao[j] = D(j)
  }
  j = j + 1
}
critic = min(tao)         #临界值
plot(ts(aipx,start = c(2001,1),end = c(2019,12),frequency = 12),ylab = '',type = 'l',
          main = '')
abline(h = 2.116788)
k1 = rep(0,k)                #变点位置
for(i in 1:k){
  j = 1
  if(D(j) == -2.116788){
    k1[j] = j
  }
  j = j + 1
}