可以认为这是R与空间计量经济学的第一部分的内容. 在上篇中介绍了空间计量经济学的一些可视化层次的内容,不过对于研究空间对象,这也是必须的。在这个第二部分,一起探究些关乎空间统计的简单内容。

1. 空间权重矩阵

1.1 以相邻定义邻居

setwd("F:/Chankey homepage work directory/spatial econometric")  
require(maptools)  
require(spdep)  
spdata <- readShapePoly("china31/region.shp")  
plot(spdata)  

nb1 <- poly2nb(spdata)  
nb1listw <- nb2listw(nb1, style = "W", zero.policy = NULL)  
nb1mat <- nb2mat(nb1, style = "W", zero.policy = NULL)  
plot(nb1, coordinates(spdata), add = TRUE, pch = ".", col = 2)  

1.2 k最近距离定义邻居

knn <- knearneigh(coordinates(spdata), k = 2, longlat = TRUE)  
plot(spdata)  
nb2 <- knn2nb(knn)  
nb2listw <- nb2listw(nb2, style = "W", zero.policy = NULL)  
nb2mat <- nb2mat(nb2, style = "W", zero.policy = NULL)  
plot(nb2, coordinates(spdata), add = TRUE, pch = ".", col = 2)  

1.3 以距离定义的邻居

knn.1n <- knearneigh(coordinates(spdata), longlat = TRUE)  
nb3 <- knn2nb(knn.1n)  
distance <- nbdists(nb3, coordinates(spdata), longlat=TRUE)  
distvec <- unlist(distance)  
mindist <- min(distvec)  
maxdist <- max(distvec)  
nb4 <- dnearneigh(coordinates(spdata), mindist, maxdist, longlat = TRUE)  
nb4listw <- nb2listw(nb4, style = "W", zero.policy = NULL)  
nb4mat <- nb2mat(nb4, style = "W", zero.policy = NULL)  
plot(spdata)  
plot(nb4, coordinates(spdata), add = TRUE, pch = ".", col = 2)  

2. Moran’s Index Computing

R的立身之本是作图和统计分析,地图自然不例外.

1. 画图

在plot函数,内部是支持地图制作的,只不过是单纯的把经纬度当作笛卡尔坐标系来绘图,而忽略了地球的其实是个球。其距离的基础是直线距离,所以plot画出来的地图相当于,球面上的中国从外向里的投影。

无论是plot还是ggmap画地图都是基于边界数据文件.shp,这个中国地理信息系统是提供的,不过貌似最近在更新,也就关闭下载链接了。不过,文尾我会提供我收藏的两套数据,中国31个省和全国(包括南海和台湾)的边界数据(后者还包括地级市乃至更详细的边界以及河流山川数据)。下文以31个省的中国地图为例。

首先, 将下载的china31文件放到R工作目录(里面必须包含.dbf,.shp,.shx,.prj三个文件,才能解析其中的.shp文件) 然后执行代码:

setwd('H:/work/spatial econometirc/lxy/china31/china31') 
 
library(maptools, quietly = TRUE) # 如果预先没有下载,请先下载   

mydat <- readShapePoly("region.shp") #读写边界文件

data <-slot(mydat, "data") # 编译出myat的内容,为data.frame

head(data)

  	SCORE    POLY_ID POPNAME     CODE      X      Y       ID   
0  2.4909     1     浙江 	330000000 120.10 29.10 Zhejiang   
1  0.5620     2     云南 	530000000 101.30 24.14   Yunnan   
2 -0.0868     3     新疆 	650000000  85.66 42.00 Xinjiang   
3 -0.7547     4     西藏 	540000000  89.12 31.10   Xizang   
4  1.2151     5     四川 	510000000 102.90 30.28  Sichuan   
5  0.3292     6     陕西 	610000000 108.76 34.12  Shaanxi   

plot(mydat)   

map6

一幅清新,简洁的中国地图就呈现了.

当然,好玩的有很多,比如merge文件提供了中国县级行政区域代码和经纬度,那可以反应到地图上来. 当你获知某行业的公司所在经纬度,那你可以在地图上反映行业集聚信息.

match <- read.csv("match.csv", header = T)
match <- match[sample(1:3000,1500),]
head(match)

		id    lon    lat   
268  210321  41.41 122.44   
1303 441502  22.78 115.37   
1058 411081  34.14 113.49  
805  360430  29.90 116.55   
2475 411224 111.05  34.05   
2643 450804 109.45  23.13   
plot(mydat) points(match[,3],match[,2], pch = 20, col = 4, cex = 0.4)

map7

另一个,关心的问题是比如全国的各省的Gini index 的水平,要将其在地图上呈现,就涉及到属性值分布图,这里提供一个函数plot.map(来自网络,未知原作者,如果冒犯,请联系我)

library(RColorBrewer) #同上,请先安装包
library(classInt)
library(maptools)

plot.map <- function(data, poly, color = NULL, ncl = 5, main = NULL){
  if(is.null(color))color<-"PuBu"
  int <- classIntervals(data, n=ncl, style="quantile")$brks
  pal <- brewer.pal(ncl, color)
  cols <- pal[findInterval(data, int, rightmost.closed=T)]
  plot(poly, col=cols)
  title(main = main)
  int <- as.character(round(int, 1))
  inttext <- vector()
  for(i in 1:(length(int)-1)){
    temp <- paste("[", int[i], "-", int[i+1], "]")
    inttext[i] <- temp
  }
  legend(130,30, fill=pal, legend=inttext, cex=0.5, title=paste("quantile:", "score"))
}

plot.map中data指各省的属性值,poly指通过readShapePoly读取的R能识别的边界数据,color值属性值的颜色(默认为PuBu色), ncl指将各省属性值分为几个层次(默认为5个层次),main为图形标题。

plot.map(1:31, mydat, ncl = 4, main = “Just for fun”)

map8

2. 空间统计分析

2.1 距离计算,geocode和revgeocode

自己辛苦画地图多累啊,何不利用现成的。现成的,优良的,还能拿来用的,当然是google map啦。ggmap包就是利用有限开源google map API来做这些事情的,并将绘图进一步用ggplot呈现,解决了上文提到的,向内投影的问题。这里不谈用ggmap来画地图,清晰爽利和拓展性不及上面详细说过的plot函数。这里说点plot不能做的。

在设计到地理经济学或者空间计量经济学等,距离是一个很重要的指标,反应了运输成本啊什么的;好吧,其实不用架这么高,如果我想知道开车从广州暨南大学去杭州要多久,发现百度不顶用后,自己算下呗。

library(ggmap) # 需要先装ggplot2 
mapdist(from = "Guangzhou Jinan university",
 		to = "Hangzhou", 
		mode = c("driving", "walking", "bicycling"))

mode指的是开车,步行和自行车的三者衡量的距离参数。mapdist函数就会调用google map里收集的数据帮你测算距离。其中,如果设计到中转站,就要计算最短距离了,当把from和to的参数改为data.frame形式的经纬度时,就可以反应出距离矩阵,然后得到最短距离。(需要提醒的是,google map API调用每日是有限制的,当然你也可以改动mapdist内部的函数,利用完全开源的其它地图的API)

如果单纯计算球面距离的话,问题就简单得多了,geosphere包可以做这个事情,而且提供了多种球面距离计算公式(不同适用对象,误差不同).

library(geosphere)
distm(x, y, fun = distHaversine)

x和y为分别为两点的经纬度(可以是data.frame), fun表示的是不同的球面距离,其中差别可以维基百科great circle distance.

回到ggmap,它还提供地址解码和反解码的功能。

比如你知道地理名称,比如Guangzhou,想知道其经纬度可以用到geocode函数(当然,并不是所有的地理位置字符,google map都能识别). 反解码,顾名思义,就是指导经纬度,想要获取其地理位置信息,用revgeocode函数。

此乃第一部分到此为止,第二部分介绍一些基础的空间计量经济的内容,比如几类空间邻居的定义,转换,空间权重矩阵的生成,moran和Geary index的计算等.

最近科学上网,似乎越来越困难了,感觉有点窒息,感觉用google能比用Baidu带来的知识积累量要快几个数量级。

其实,如果你google一下,在R里如何实现并行计算,或者搜索foreach包。帮助文档和相关介绍博客何其多,精品也实在不少,我也是通过这些内容和不断google之,才学到这点知识的。那既然如此,这篇blog还有存在的必要吗?有,必然有,可能意义更多的集中在于我自身。从思想上让我写博客的是刘未鹏老师的两篇文章为什么你应该(从现在开始就)写博客书写是为了更好的思考.技术问题,是通过陈偃平老师的Github Pages极简教程解决的.水到渠成……

本文的运行和技术环境是windows 32位,两核操作系统.

R在前几年,被人诟病的最多的就是计算性能问题。不过随着这几年,像ff,bigmemory,doparallel和dosnow等相关包的出现,算是部分解决了直接制约R发展的几个问题,前者是大数据存储和操作,后者是并行计算。

1. 什么是并行运算

概念介绍,首先wiki,这里谈及的是空间多核并行。

并行运算是相对串行运算而言的,正如初中学过的并行和串行电路。前者指的是一次执行多个算法或任务,以提高运算速度。这里的执行多个算法就是基于多个处理器并发计算。所以具体同时实行几个任务,取决于你电脑的内核数(而不是32位还是64位-_-!)。或许你会问,我怎么知道我电脑几核的啊?请ctrl+alt+delete打开任务管理器,点击性能栏,出现几个cpu使用记录,就是几核了,一般笔记本为两核无误。
5

2. 怎么添加并行运算

我在zk的论文项目里,设计的并行运算主要用到foreach和doParallel包,这是在win系统,在linux系统里后者是doMc包。这是由旨在提高R操作性能的公司R Evolution Analytics贡献给R社区的。

Getting Started with doParallel and foreach这篇文档,应该是直接的入门材料,由R Evolution Analytics公司的 Steve WestonRich Calaway 提供。

library(foreach)
library(doParallel)
cl <- makeCluster(2) # 2为你电脑的核数
registerDoParallel(cl)

foreach(i = 1:10) %dopar% {  
	loop body
}

stopCluster(cl) # 运行结束后可关闭分支

foreach默认返回的结果为list,每一个循环保存为list的一个元素。当然,该包提供了很多其它的形式,甚至是自定义的形式,比如:

foreach(i = 1:10,  .combine=cbind) %dopar% {  
	loop body
}

就是列并的形式。这在实际问题处理中,提供了巨大的方便。

3. 并行运算有哪些不好的地方?

怎么样,打开任务管理器看看性能栏,是不是cpu爆满利用率100%,看着它满负荷的为你当牛做马,心里是不是爽歪歪!

需要指出的是,并行虽然可以提高运行速度,但是

  • cpu满载,对电脑总是没那么好(四核,八核的实验室台式电脑才没这么心疼呢)
  • 不能像for循环可以中断保存信息
  • 通过message或者cat函数设计进度条,无效

另外,在我们的plyr包中,ply类的函数里都提供了并行形式,如虎添翼啊!
简单的并行introduction到此为止,下面不妨来看看并行提高运行速度的实际效果:

require(foreach, quietly = TRUE)  
require(doParallel, quietly = TRUE)  
cl <- makeCluster(2)  
registerDoParallel(cl)  

t0 <- Sys.time()  
list1 <- foreach(i = 1:3000) %do% {  
  rn <- rnorm(30000)  
  sample(rn, 10*i)  
}  
t1 <- Sys.time()  
t1 - t0  

t0 <- Sys.time()  
list1 <- foreach(i = 1:3000) %dopar% {  
  rn <- rnorm(30000)   
  sample(rn, 10*i)  
}  
t1 <- Sys.time()  
t1 - t0   

比较两个t1-t0的运行结果发现,也没快几秒啊。逗我呢?必须指出的是,并行运算在计算小程序时,并不会显著的提高速度,有时候反而会慢几秒。这有种杀鸡焉用牛刀的感觉,不过诚然,在较大程序时,至少10分钟以上的程序,并行的优良性还是显著的。(我会在我的破本本上做实验?分分钟自动关机给你看……)

上个月处理一些医学统计的冗杂工作,提供的都是30多个工作簿,而且每个工作簿里都有30多个工作表。需要对其进行一些简单的统计分析,用R处理要方便很多。自然,问题就是多sheets excel在R中的输入和输出问题。

1. 用VBA函数解决多excel合并的问题

我们需要处理的文件夹里有30多个excel文件,因为它们的同质性,尝试将其以sheets的形式合并在一个excel中。

合并工作簿

Sub 工作薄间工作表合并()  
Dim FileOpen  
Dim X As Integer  
Application.ScreenUpdating = False  
FileOpen = Application.GetOpenFilename(FileFilter:="Microsoft Excel文件(*.xlsx),*.xlsx",MultiSelect:=True, Title:="合并工作薄")  
X = 1   
While X <= UBound(FileOpen)    
Workbooks.Open Filename:=FileOpen(X)   
Sheets().Move After:=ThisWorkbook.Sheets(ThisWorkbook.Sheets.Count)   
X = X + 1   
Wend    
ExitHandler:    
Application.ScreenUpdating = True    
Exit Sub   

errhadler:   
MsgBox Err.Description   
End Sub   

在要合并的30多个工作簿的文件夹中,新建一个excel,打开后,在sheet1的位置,右击查看代码,然后黏贴运行。提示你打开需要合并的excel工作簿,将其全选点击打开即可。

2. 用VBA函数解决多sheets合并为一个sheet的问题

在将excel合并之前或者合并之后,excel中若含有多数sheets,将其合并到同一个sheet的VBA代码:

Sub 合并当前工作簿下的所有工作表()
Application.ScreenUpdating = False
For j = 1 To Sheets.Count
   		If Sheets(j).Name <> ActiveSheet.Name Then
   		X = Range("A65536").End(xlUp).Row + 1
   		Sheets(j).UsedRange.Copy Cells(X, 1)
   		End If
Next
Range("B1").Select
Application.ScreenUpdating = True
MsgBox "当前工作簿下的全部工作表已经合并完毕!", vbInformation, "提示"
End Sub

这里需要注意的是要确保sheets的格式要一样。比如很自然的一件事是,在合并之前,从第二个sheet起,就无表头标题存在的必要。一个简单的小方法是,在sheet2处右击,选定全部工作表之后,对某一sheet操作就可以了。也可以在选定所有工作表后,按住ctrl后点击sheet1,这样操作就排除了sheet1,可以保留其表头和格式。

3. 用自定义函数读出list对象到多sheets的表格

在office里,我唯一钟情的只有excel,因其及其灵活、快速的单元格操作,加之VBA,千变万化,及其便易。经常碰到在R里面难以处理的问题,或者实现稍微复杂的问题,不妨通过excel和VBA去解决。

普通的数据读出,write函数或者xlsx包里的wirte.xlsx的函数都没问题。但是在R里在处理同质性问题时,经常将其处理为list和array对象。下面提供一个函数,可以将list对象读出为多sheets的工作簿:

require(xlsx, quietly = TRUE)

save.xlsx <- function (file, objects) {
  		nobjects <- length(objects)
  		for (i in 1:nobjects) {
  		  if (i == 1)
  		write.xlsx(objects[[i]], file, sheetName = names(objects)[i])
	  else write.xlsx(objects[[i]], file, sheetName = names(objects)[i], append = TRUE)
  		}
	print(paste("Workbook", file, "has", nobjects, "worksheets."))
}

save.xlsx('file.xlsx', object) 

其中file为读出的excel名称,可以定义输出位置,object为你需要输出的list对象。sheet的名称就是list元素的名称。
通过xlsx包,貌似是需要调用javascript的,当list的对象较大时,会经常报错“error in:js out of memory”。我也不清楚具体的原因,不过我们可以采用分而治之的办法。

names <- c()
for (i in 1:18) { # 划分为18个excel
  	  names <- c(names, paste('list', i,'.', 'xlsx', sep = ''))  # excel的名字取为list1.xlsx,list2.xlsx类推
}

for ( i in 1:18) {
  save.xlsx(names[i], companydata.dest.list[((10*i)-9):(10*i)]) # 以list的每10个元素读入为一个excel,共180个元素
}

1.状态空间模型

状态空间模型(state-space-model),又称作隐马尔科夫模型(Hide Markov model).它描述的是一类广泛、灵活的时间序列模型。时间序列里诸多经典模型(ARMA,ARIMA,GARCH)等皆为状态空间模型的特殊情况。其试图以系统的观点来看待研究对象和其所处环境。而用来描述系统的概念——状态和信号。
状态序列:$S_{0:T} ={S_0, S_1, \cdots, S_T }$

可观察序列: $Y_{0:T} = {Y_0, Y_1, \cdots, Y_T }$

状态是无法直接观察到的,它是系统的真实描述;我们只能通过信号来部分推测系统的真实状态。故状态空间模型分为两个部分:

初始状态:

进化过程:

可观察过程:

其中$\theta$为模型的参数.

状态空间模型应用及其广泛,比如目标定位以及动态追踪问题,信号传输的降噪问题等等;在经济金融领域,比如对真实CPI的滤波求解,在高频以及超高频时间序列中真实价格序列的求解问题,这里拿经典的随机波动率模型为例,$X_1 \sim N(0,\frac{\sigma^2}{1-\alpha^2})$

其中,通过error项的正态性,可以得出.

那么问题来了,已知, 依赖此模型到底能得到哪些感兴趣的东西?

  • 滤波问题:对进行推断,以及从含噪的观察序列中怎么得出真实的状态序列;过滤噪音,滤波名称由此而来.(动态滤波)
  • 预测问题:对进行推断,在现在掌握的所有信息下,对下个状态进行预测.(目标跟踪)
  • 平滑问题:对进行推断,当获知全部可观察信息时,能不能更好的对状态变量进行估计.(信号传输)
  • 参数学习:对进行推断,怎么根据已知观察序列信息,对模型参数进行估计.

其实状态空间模型,只是提供了一个灵活的框架。根据你的问题,设计模型,然后对上述四个问题进行解答。
在滤波问题上,最经典的乃Kalman滤波,在状态空间模型为线性和高斯的情形下,能得到精确解,这一点是伟大的。然而,非线性才是自然的本质,在原来Kalman滤波的框架下,虽然提出了UKF和EKF(通过高斯近似和线性近似的方法)解决非线性和非高斯的问题,可是效果不太好。最新几年,粒子滤波被提出来解决这个问题(这是后续blog的重点)。

2.微观结构噪音

2.1 低频,高频和超高频数据

低频时间序列数据,是通常我们见到和研究的数据,其抽样频率一般为天,周,月等,比如每天的收盘价. 而高频时间序列数据虽然也是抽样数据,但是其抽样频率更低,聚焦在日内(秒,分,小时). 而超高频时间序列数据,就完全是根据交易实际发生所记录的数据(tick by tick),国外顶尖团队已达微妙级,国内为毫秒级. 这类数据很明显有两类特征:

  • 交易时间的随机性

  • 交易价格被微观结构噪音严重污染

2.2 微观结构噪音描述

微观结构噪音,主要包括以下三个方面:

  • 离散性: 最小变动单位, 由交易所设定的交易规则 $\frac{1}{8},\frac{1}{32},\frac{1}{100}$

  • 集聚性: Harris(1991) 发现, (以$\frac{1}{8}$ 为例),

    整数 > 小数部分为0.5 > 小数部分为0.25或0.75 > 小数部分为$\frac{1}{8}$的奇数倍

    并通过计量工具证明此关系的显著性,并将原因解释为降低交易协商成本

    cdd

  • 非集聚性: _Black(1986)_指出:固定交易成本, 做市商的库存控制, 非对称信息交易;这解释了离异值的发生

2.3 微观结构噪音的影响

  • Bandi & Russell (2006): 已实现波动率的过高估计
  • Asparouhona (2010): 普通最小二乘回归下,得出有偏且不一致的系数估计量
  • Duan & Fulop (2009): 在信用风险结构模型中,对企业资产波动率显著的过高估计

结论:
Hasbrouck (1996)指出市场信息对证券价格起到永久性的影响而噪音仅起到短暂性影响. 由于微观结构噪音的存在,传统针对低频时间序列数据的建模方法不适用于超高频时间数据

2.4 微观结构噪音建模

对上述微观结构噪音分别建模为:

  • 离散噪音: $Round[S(t_i),\frac{1}{M}]$,$\frac{1}{M}$为最小变动单位.

  • 非集聚噪音: 由服从双几何分布的随机变量刻画的随机游走 $Y^{‘}(t_i)=Round[S(t_i),\frac{1}{M}] +\nu_i$

  • 集聚型噪音,(设定 M = 8 为例):

    • 如果 $Y^{‘}(t_i)$ 分数部分是 $\frac{1}{8}$ 的偶数倍, $Y(t_i)= Y^{‘}(t_i)$

    • 否则,以概率$\alpha$偏移$Y^{‘}(t_i)$ 到最近$\frac{1}{4}$的奇数倍

      以概率$\beta$偏移$Y^{‘}(t_i)$ 到最近$\frac{1}{4}$的偶数倍

      以概率$\gamma$偏移$Y^{‘}(t_i)$ 到最近整数

      以概率$1-\alpha - \beta - \gamma$保留$Y(t_i) = Y^{‘}(t_i)$

    参数 $\alpha,\beta,\gamma$ 可以由历史数据通过频率的方法求出,Zeng(2003).

设定$\rho=0.4$, 由实际超高频价格序列的集聚性和由该模型模拟的集聚性对比
dcd
—————————————————————–

在超高频时间序列数据中,真实价格序列是被微观结构噪音严重污染的,故而通过状态空间模型进行处理。这里花大部分笔墨刻画微观结构噪音,其实描述的状态空间模型的第二部分,即可观察过程. 后续文章,进一步描述状态空间模型在研究超高频时间序列数据问题的第一部分——状态进化过程.