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

2024-04-22 10:38
文章标签 算法 语言 程序 方法 老师 计量 经济 分析方法 参照 景气 张永军

本文主要是介绍经济景气计量分析方法（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  #固有模式函数
    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)
}

这篇关于经济景气计量分析方法（R语言程序算法）参照张永军老师的方法的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/925561。

相关文章

Nginx设置连接超时并进行测试的方法步骤

Nginx设置连接超时并进行测试的方法步骤

《Nginx设置连接超时并进行测试的方法步骤》在高并发场景下,如果客户端与服务器的连接长时间未响应,会占用大量的系统资源,影响其他正常请求的处理效率,为了解决这个问题,可以通过设置Nginx的连接... 目录设置连接超时目的操作步骤测试连接超时测试方法：总结：设置连接超时目的设置客户端与服务器之间的连接
阅读更多...
Java判断多个时间段是否重合的方法小结

Java判断多个时间段是否重合的方法小结

《Java判断多个时间段是否重合的方法小结》这篇文章主要为大家详细介绍了Java中判断多个时间段是否重合的方法,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 目录判断多个时间段是否有间隔判断时间段集合是否与某时间段重合判断多个时间段是否有间隔实体类内容public class D
阅读更多...
Python使用国内镜像加速pip安装的方法讲解

Python使用国内镜像加速pip安装的方法讲解

《Python使用国内镜像加速pip安装的方法讲解》在Python开发中,pip是一个非常重要的工具,用于安装和管理Python的第三方库,然而,在国内使用pip安装依赖时,往往会因为网络问题而导致速... 目录一、pip 工具简介1. 什么是 pip？2. 什么是 -i 参数？二、国内镜像源的选择三、如何
阅读更多...
IDEA编译报错“java: 常量字符串过长”的原因及解决方法

IDEA编译报错“java: 常量字符串过长”的原因及解决方法

《IDEA编译报错“java:常量字符串过长”的原因及解决方法》今天在开发过程中,由于尝试将一个文件的Base64字符串设置为常量,结果导致IDEA编译的时候出现了如下报错java:常量字符串过长,... 目录一、问题描述二、问题原因2.1 理论角度2.2 源码角度三、解决方案解决方案①：StringBui
阅读更多...
Linux使用nload监控网络流量的方法

Linux使用nload监控网络流量的方法

《Linux使用nload监控网络流量的方法》Linux中的nload命令是一个用于实时监控网络流量的工具,它提供了传入和传出流量的可视化表示,帮助用户一目了然地了解网络活动,本文给大家介绍了Linu... 目录简介安装示例用法基础用法指定网络接口限制显示特定流量类型指定刷新率设置流量速率的显示单位监控多个
阅读更多...
Java覆盖第三方jar包中的某一个类的实现方法

Java覆盖第三方jar包中的某一个类的实现方法

《Java覆盖第三方jar包中的某一个类的实现方法》在我们日常的开发中,经常需要使用第三方的jar包,有时候我们会发现第三方的jar包中的某一个类有问题,或者我们需要定制化修改其中的逻辑,那么应该如何... 目录一、需求描述二、示例描述三、操作步骤四、验证结果五、实现原理一、需求描述需求描述如下：需要在
阅读更多...
JavaScript中的reduce方法执行过程、使用场景及进阶用法

JavaScript中的reduce方法执行过程、使用场景及进阶用法

《JavaScript中的reduce方法执行过程、使用场景及进阶用法》：本文主要介绍JavaScript中的reduce方法执行过程、使用场景及进阶用法的相关资料,reduce是JavaScri... 目录1. 什么是reduce2. reduce语法2.1 语法2.2 参数说明3. reduce执行过程
阅读更多...
C#中读取XML文件的四种常用方法

C#中读取XML文件的四种常用方法

《C#中读取XML文件的四种常用方法》Xml是Internet环境中跨平台的,依赖于内容的技术,是当前处理结构化文档信息的有力工具,下面我们就来看看C#中读取XML文件的方法都有哪些吧... 目录XML简介格式C#读取XML文件方法使用XmlDocument使用XmlTextReader/XmlTextWr
阅读更多...
python使用fastapi实现多语言国际化的操作指南

python使用fastapi实现多语言国际化的操作指南

《python使用fastapi实现多语言国际化的操作指南》本文介绍了使用Python和FastAPI实现多语言国际化的操作指南,包括多语言架构技术栈、翻译管理、前端本地化、语言切换机制以及常见陷阱和... 目录多语言国际化实现指南项目多语言架构技术栈目录结构翻译工作流1. 翻译数据存储2. 翻译生成脚本
阅读更多...
C++初始化数组的几种常见方法(简单易懂)

C++初始化数组的几种常见方法(简单易懂)

《C++初始化数组的几种常见方法(简单易懂)》本文介绍了C++中数组的初始化方法,包括一维数组和二维数组的初始化,以及用new动态初始化数组,在C++11及以上版本中,还提供了使用std::array... 目录1、初始化一维数组1.1、使用列表初始化（推荐方式）1.2、初始化部分列表1.3、使用std::
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2