经济景气计量分析方法（R语言程序算法）参照张永军老师的方法

#==========方 法==========#

#1.经济指标时间序列数据的分解

#2.先行、一致、滞后指标的选取

#3.景气指数的计算

#4.景气指数的检验与修正（该步未写相关算法）
#=========================#

#首先肯定是读取数据，假设读进来之后命名为data1，而且数据的第一列是序号或者日期

data2 <- na.omit(data1[,-1]) #第一列是序号

n <- nrow(data2) #样本量

#1.经济指标时间序列数据的分解（采用的方法是stl分解）

#1.1季节调整

period<-11 #周期

fstl <- function(x){

  index1.ts <- ts(as.double(x),frequency = 11,start = 2009+2/12)#构造'ts'结构时间序列

  index1.ts.stl1 <- stl(index1.ts,t.window=13,s.window=ceiling((1.5*period)/(1-(1.5/13))),robust = T)

  return(index1.ts.stl1)

} 

stl1 <- apply(data2,2,fstl)

data3 <- c()

for(i in 1:length(stl1))

  data3 <- cbind(data3, stl1[[i]]$time.series[,2])

colnames(data3) <- colnames(data2)

#2.先行、一致、滞后指标

#KL信息量（有的称为KL距离或者KL信息熵） 

fmm <- function(x) { #数据归一化

  (x - min(x))/(max(x) - min(x)) + 1 #加1是为了不出现0

}

datamm <- apply(data3, 2, fmm) #归一化，保证每个值都大于0

fp <- function(x) { #得到概率分布

  x/sum(x)

}

datap <- apply(datamm, 2, fp) #概率分布

fkl <- function(p1,p2,L=4) { #L是最大滞后期

  k <- c()

  for(l in -L:L) {

    k1 <- 0

    if(l > 0) {

      for(j in 1:(length(p1)-l))

        k1 <- k1 + p1[j]*log(p1[j]/p2[j+l])

    } else {

      for(j in (abs(l)+1):length(p1))

        k1 <- k1 + p1[j]*log(p1[j]/p2[j+l])

    }

    k <-c(k, k1) 

  }

  c(-L:L)[which.min(abs(k))] #KL信息量接近于0时对应的滞后

}

lag1 <- matrix(0, ncol(datap), ncol(datap)) #变量间的先行、一致、滞后关系,保证对角线为0

for(i in 1:ncol(datap)) {

  for(j in c(1:ncol(datap))[-i])

  lag1[j, i] <- fkl(datap[,i],datap[,j])

}

#EMD和符号化相结合（该方法是借鉴文献：基于EMD与K_means算法的时间序列聚类_刘慧婷）

#方法说明：

#1.先用EMD分解,得到数据的趋势序列

#2.计算趋势序列的符号化距离

library(EMD)

#library(hht)

library(TSclust)

fEMDSAx <- function(var1, var2, L=4) {

  fEMD <-function(var1, var2){

    #EMD分解,得到数据的趋势序列

    n <- length(var1)

    res_all <- matrix(NA,2,n)

    em1 <- emd(var1, tol = sd(var1)*0.01^2, boundary = "wave") #符号化距离，得到emd分解

    #imf1 <- em1$imf

    res1 <- em1$residue

    res_all[1,] <- res1

    

    em2<-emd(var2, tol = sd(var2)*0.01^2, boundary = "wave")

    #imf2 <- em2$imf  #固有模式函数

4000

    res2 <- em2$residue #趋势序列

    res_all[2,] <- res2

    m <- max(em1$nimf,em2$nimf) #最大的分隔数

    

    #计算趋势序列的符号化距离

    if((m^2)<n) k <- m^2 else k <- sqrt(n)

    d2<-diss(res_all, "MINDIST.SAX", k, alpha = k) #dist类

    d2

    

  }

  

  for(l in -L:L) {

    k1 <- 0

    if(l > 0) {

      k1 <- fEMD(var1[1:(length(var1)-l)], var2[(1:(length(var1)-l))+l])

    } else { 

      k1 <- fEMD(var1[(abs(l)+1):length(var1)], var2[((abs(l)+1):length(var1))+l])

    }

    k <-c(k, k1) 

  }

  c(-L:L)[which.min(abs(k))] 

}

lag2 <- matrix(0, ncol(datamm), ncol(datamm)) #变量间的先行、一致、滞后关系,保证对角线为0

for(i in 1:ncol(datamm)) {

  for(j in c(1:ncol(datamm))[-i])

    lag2[j, i] <- fkl(datamm[,i],datamm[,j])

}

#3.权重的确定

#3.1.熵值法

fentropy <- function(x) {

  n4 <- dim(x)

  ej1 <- c()

  for(j in 1:n4[2]) {

    a1 <- 0

    for(i in 1:n4[1]) {

      a2 <- x[i, j]*log(x[i, j])

      a1 <- a1 + a2

    }

    ej1 <- c(ej1, a1)

  }

  ej <- c()

  for(i in 1:n4[2]) {

    ej2 <- (-1)/(nrow(x)*ej1[i])

    ej <- c(ej, ej2)

  }

  w <- (1-ej)/sum(1-ej)

  w

}

#3.2.层次分析法(AHP)

#先建立判断矩阵，然后调用ahp函数得到权重

library(stats4)

library(pmr)

fAHP <- function(data) {

  n <- ncol(data)

  A <- matrix(1, n, n) #判断矩阵

  for(i in 1:(n-1))

    for(j in (i+1):n) {

      A[j, i] <- abs(cor(data[,i], data[,j]))

      A[i, j] <-  1/A[j, i] 

    }

  ahp(A) #A是判断矩阵

}

ahp1 <- fAHP(datamm)

#4.扩散指数DI

fDI <- function(data, j=3) { #j为基期

  N <- ncol(data)

  w <- as.numeric(fentropy(data))

  DI <- 0

  for(i in 1:N){

    bool1 <- 0

    for(k in (j+1):nrow(data)){

      bool1 <- bool1 + (data[k,i]>data[k-j,i])

    }

    num <- 0.5*ceiling(0.5*nrow(data))

    if(bool1>num) I <- 1 else if(bool1==num) I <- 0.5 else I <- 0 

    DI <- DI + w[i]*I

  }

  DI

}

di1 <- fDI(datamm);di1

#5.合成指数CI

fCI <- function(lag1, vn=1) { #vn:所要研究变量（因变量）所在的列号

  p1 <- as.matrix(datamm[,which(lag1[,vn]<0)])#先行指标组

  p2 <- as.matrix(datamm[,which(lag1[,vn]==0)])#一致指标组

  p3 <- as.matrix(datamm[,which(lag1[,vn]>0)])#滞后指标组

  if(length(p1)>0){ #存在先行指标

    n11 <- dim(p1) 

    c1 <-  array(0,n11)

    for(i in 1:(n11[1]-1))

      for(j in 1:n11[2])

        c1[i, j] <- p1[i+1, j] - p1[i, j]

  }

  if(length(p2)>0) { #存在一致指标

    n22 <- dim(p2)

    c2 <-  array(0,n22)

    for(i in 1:(n22[1]-1))

      for(j in 1:n22[2])

        c2[i, j] <- p2[i+1, j] - p2[i, j]

    

  }

  if(length(p3)>0) { #存在滞后指标

    n33 <- dim(p3)

    c3 <-  array(0,n33)

    for(i in 1:(n33[1]-1))

      for(j in 1:n33[2])

        c3[i, j] <- p3[i+1, j] - p3[i, j]

  }

  if(length(p1)>0&&length(p2)>0&&length(p3)>0) c123 <- cbind(c1, c2, c3) else 

    if(length(p1)<=0&&length(p2)>0&&length(p3)>0) c123 <- cbind(c2, c3) else 

      if(length(p1)>0&&length(p2)<=0&&length(p3)>0) c123 <- cbind(c1, c3) else c123 <- cbind(c1, c2)

  a <- matrix(0,1,ncol(c123))

  for(i in 1:ncol(a))

    for(j in 1:ncol(c123))

      a[i] <- sum(abs(c123[, j]))/(n-1)

  s <- array(0,dim(c123))

  for(i in 1:ncol(c123))

    s[,i] <- c123[, i]/a[i]

  if(length(p1)>0&&length(p2)>0&&length(p3)>0) { #存在先行、一致、滞后指标

    w1 <- fentropy(p1)

    w2 <- fentropy(p2)

    w3 <- fentropy(p3)

    s1 <- as.matrix(s[,1:ncol(p1)])

    s2 <- as.matrix(s[,(ncol(p1)+1):(ncol(p1)+ncol(p2))])

    s3 <- as.matrix(s[,(ncol(p1)+ncol(p2)+1):ncol(s)])

    r1 <- r2 <- r3 <- matrix(0, 1, n)

    for(i in 1:(n-1)) {

      r1[i] <- sum(s1[i,]*w1)/sum(w1)

      r2[i] <- sum(s2[i,]*w2)/sum(w2)

      r3[i] <- sum(s3[i,]*w3)/sum(w3)

    }

    r <- rbind(r1, r2, r3)

    f <- matrix(0, 1, 3)

    for(i in 1:3)

      f[i] <- (sum(abs(r[i,]))/(n-1))/(sum(abs(r[2,]))/(n-1))

    v1 <- array(0,dim(r))

    for(i in 1:nrow(r))

      v1[i,] <- r[i,]/f[i]

    I <- matrix(0, 3, n)

    I[, 1] <- 100

    for(i in 1:3)

      for(j in 2:n)

        I[i, j] <- I[i, j-1]*(200+v1[i, j-1])/(200-v1[i, j-1])

    I <- I/100

    CI <- matrix(0, 3, n)

    for(i in 1:3) {

      I2 <- sum(I[i,])/n

      CI[i, ] <- (I[i,]/I2)*100

      if(i==1) plot(CI[i,],type = "l",col=1) else lines(CI[i,],col=i)

    }

    legend("topleft", cex=0.7, col = c(1:3), lty = 1, legend = c("先行","一致","滞后"))

  } else if(length(p1)<=0&&length(p2)>0&&length(p3)>0) { #存在一致、滞后指标

    w2 <- fentropy(p2)

    w3 <- fentropy(p3)

    s2 <- as.matrix(s[,1:ncol(p2)])

    s3 <- as.matrix(s[,(ncol(p2)+1):ncol(s)])

    r2 <- r3 <- matrix(0, 1, n)

    for(i in 1:(n-1)) {

      r2[i] <- sum(s2[i,]*w2)/sum(w2)

      r3[i] <- sum(s3[i,]*w3)/sum(w3)

    }

    r <- rbind(r2, r3)

    f <- matrix(0, 1, 2)

    for(i in 1:2)

      f[i] <- (sum(abs(r[i,]))/(n-1))/(sum(abs(r[2,]))/(n-1))

    v1 <- array(0,dim(r))

    for(i in 1:nrow(r))

      v1[i,] <- r[i,]/f[i]

    I <- matrix(0, 2, n)

    I[, 1] <- 100

    for(i in 1:2)

      for(j in 2:n)

        I[i, j] <- I[i, j-1]*(200+v1[i, j-1])/(200-v1[i, j-1])

    I <- I/100

    CI <- matrix(0, 2, n)

    for(i in 1:2) {

      I2 <- sum(I[i,])/n

      CI[i, ] <- (I[i,]/I2)*100

      if(i==1) plot(CI[i,],type = "l",col=1) else lines(CI[i,],col=i)

    }

    legend("topleft", cex=0.7, col = c(1:2), lty = 1, legend = c("一致","滞后"))

  } else if(length(p1)>0&&length(p2)<=0&&length(p3)>0) { #存在先行、滞后指标

    w1 <- fentropy(p1)

    w3 <- fentropy(p3)

    s1 <- as.matrix(s[,1:ncol(p1)])

    s3 <- as.matrix(s[,(ncol(p1)+1):ncol(s)])

    r1 <- r3 <- matrix(0, 1, n)

    for(i in 1:(n-1)) {

      r1[i] <- sum(s1[i,]*w1)/sum(w1)

      r3[i] <- sum(s3[i,]*w3)/sum(w3)

    }

    r <- rbind(r1, r3)

    f <- matrix(0, 1, 2)

    for(i in 1:2)

      f[i] <- (sum(abs(r[i,]))/(n-1))/(sum(abs(r[2,]))/(n-1))

    v1 <- array(0,dim(r))

    for(i in 1:nrow(r))

      v1[i,] <- r[i,]/f[i]

    I <- matrix(0, 2, n)

    I[, 1] <- 100

    for(i in 1:2)

      for(j in 2:n)

        I[i, j] <- I[i, j-1]*(200+v1[i, j-1])/(200-v1[i, j-1])

    I <- I/100

    CI <- matrix(0, 2, n)

    for(i in 1:2) {

      I2 <- sum(I[i,])/n

      CI[i, ] <- (I[i,]/I2)*100

      if(i==1) plot(CI[i,],type = "l",col=1) else lines(CI[i,],col=i)

    }

    legend("topleft", cex=0.7,col = c(1:2), lty = 1, legend = c("先行", "滞后"))

  } else { #存在先行、一致指标

    w1 <- fentropy(p1)

    w2 <- fentropy(p2)

    s1 <- as.matrix(s[,1:ncol(p1)])

    s2 <- as.matrix(s[,(ncol(p1)+1):ncol(s)])

    r1 <- r2 <- matrix(0, 1, n)

    for(i in 1:(n-1)) {

      r1[i] <- sum(s1[i,]*w1)/sum(w1)

      r2[i] <- sum(s2[i,]*w2)/sum(w2)

    }

    r <- rbind(r1, r2)

    f <- matrix(0, 1, 2)

    for(i in 1:2)

      f[i] <- (sum(abs(r[i,]))/(n-1))/(sum(abs(r[2,]))/(n-1))

    v1 <- array(0,dim(r))

    for(i in 1:nrow(r))

      v1[i,] <- r[i,]/f[i]

    I <- matrix(0, 2, n)

    I[, 1] <- 100

    for(i in 1:2)

      for(j in 2:n)

        I[i, j] <- I[i, j-1]*(200+v1[i, j-1])/(200-v1[i, j-1])

    I <- I/100

    CI <- matrix(0, 2, n)

    for(i in 1:2) {

      I2 <- sum(I[i,])/n

      CI[i, ] <- (I[i,]/I2)*100

      if(i==1) plot(CI[i,],type = "l",col=1) else lines(CI[i,],col=i)

    }

    legend("topleft", cex=0.7, col = c(1:2), lty = 1, legend = c("先行","一致"))

  }

  CI

}

#KL信息量得到的合成指数

ci1 <- vector("list", ncol(lag1))

for(i in 1:ncol(lag1)) {

  ci1[[i]] <- fCI(lag1, vn=i)

}

#EMD_SAX得到的合成指数

ci2 <- vector("list", ncol(lag2))

for(i in 1:ncol(lag2)) {

  ci2[[i]] <- fCI(lag2, vn=i)

}