@fanxy
2018-04-23T15:52:22.000000Z
字数 3991
阅读 2731
樊潇彦
复旦大学经济学院
数量经济学
setwd("D:\\...")
rm(list=ls())
install.packages(c("nleqslv","deSolve","scatterplot3d"))
install.packages(c("zoo","quantmod","tseries","fBasics","forecast")) # 时间序列常用包
library(dplyr)
library(tidyr)
library(ggplot2)
library(nleqslv)
library(deSolve)
library(scatterplot3d)
delta=0.05; s=0.4; A=1; alpha=0.3 # 设定参数
K=1; times = seq(0, 100, by = 1) # 初始值和时间轴
# 求稳态资本 Kstar
FunKstar<-function(Kstar) {
error=s*A*Kstar^alpha - delta*Kstar # K_t+1 - K_t =0
return(error)
}
library(nleqslv) # 调用程序包,内含命令nleqslv求解非线性方程(组)
Kstar=nleqslv(5,FunKstar)$x
# 解差分方程
disK=rep(K,length(times))
for (t in 2:length(times)) disK[t]=s*A*disK[t-1]^alpha + (1-delta)*disK[t-1]
disK=data.frame(times=times,`差分方程`=disK)
# 解微分方程
Kfun=function(t,K,par){
with(as.list(c(K,par)),{
dK= s*A*K^alpha - delta*K
list(dK)
})
}
conK= ode(y=c(K=1), times=times, func=Kfun, # ode:程序包deSolve中命令,解常微分方程
parms=c(delta=delta,s=s,A=A,alpha=alpha))
# 合并数据,作图
Capital=data.frame(times=times, `微分方程`=conK[,-1])%>%
left_join(disK,by="times")
tail(Capital,5) # 最后5期结果比较
ggplot(Capital%>%gather(var,K,-times), aes(times,K, color=var))+geom_line()+
geom_hline(yintercept=Kstar,linetype = "dotdash",col="red")+
labs(title="Solow-Swan模型",x="时期",y="人均资本")+
guides(color=guide_legend(title=NULL))+theme_bw()+theme(legend.position="bottom")
parameters <- c(a = -8/3, b = -10, c = 28)
state <- c(X = 1, Y = 1, Z = 1)
times <- seq(0, 100, by = 0.01)
Lorenz <- function(t, state, parameters) {
with(as.list(c(state, parameters)), {
dX <- a * X + Y * Z
dY <- b * (Y - Z)
dZ <- -X * Y + c * Y - Z
list(c(dX, dY, dZ))
})
}
out <- ode(y = state, times = times, func = Lorenz, parms = parameters)
scatterplot3d(out[,-1], type = "l", highlight.3d=TRUE, col.axis="blue", col.grid="lightblue")
# 自已写一个可以计算矩阵的幂的函数
mat_power = function(A, n){
Apower=A
for (i in 2:n) Apower= Apower %*% A
return(Apower)
}
P=matrix(c(0.95,0.05,0.03,0.97),nrow=2,byrow=T) # 人口迁移的转移矩阵
P_star=mat_power(P,300) # 计算稳态转移矩阵
round(eigen(t(P))$values,2) # 求特征值
ev=eigen(t(P))$vectors # 求左特征向量
x=ev[,1]/sum(ev[,1]) # 将特征值为1的特征向量正规化
# 比较三种稳态人口比例
x
x%*%P
P_star[1,]
学校电子图书馆有很多中文数据库,下面以国泰安中季度国内生产总值为例,说明宏观经济数据的特征。
Qgdp=read.table("CME_Qgdp.txt", header = T, stringsAsFactors =F)
# summary(Qgdp); str(Qgdp) # 查看数据结构
library(zoo)
Qgdp=Qgdp%>%mutate(time=as.yearmon(Staper))%>%
rename(`季度GDP(累计值)`=Gdpq0101)%>%rename(`同比增长率(名义)`=Gdpq0105)%>%
select(time,`季度GDP(累计值)`,`同比增长率(名义)`)
ggplot(Qgdp%>%gather(var,value,-time),aes(time,value))+
geom_line()+facet_wrap(~var,scales="free")+
labs(title="中国宏观经济数据",x="",y="")+theme_bw()
用quantmod
包中的getSymbols()
函数从yahoo财经下载上证综指,更多指标可以查看相关网站。
library(quantmod)
ssec<-getSymbols("^SSEC",from = "1991-07-15",to = Sys.Date(),src = "yahoo",auto.assign=FALSE)
dim(ssec)
tail(ssec)
plot(ssec$SSEC.Adjusted, main = "上证综指(SSEC)")
ssec.m<-to.monthly(ssec) # 转换成月度数据
chartSeries(ssec.m,type="candlesticks",theme="white", name = "上证综指月线图",time.scale = "",up.col="red",dn.col="green")
ln_ssec=log(ssec.m$ssec.Adjusted)
# 计算月收益率,由于月收益率可能较大,不用diff(log(x))
ssec.m$ret=diff(ssec.m$ssec.Adjusted) / lag(ssec.m$ssec.Adjusted) * 100
ret=ssec.m$ret[-1,] # 去掉第一行缺失值
plot(ret, main="上证综指月收益率(%)")
library(tseries)
pp.test(ln_ssec) # DF检验,H0:有单位根
pp.test(ret)
adf.test(ln_ssec) # ADF检验,H0:有单位根
adf.test(ret)
kpss.test(ln_ssec) # KPSS检验,H0:平稳
kpss.test(ret)
实际应用中,forecast
包中的 auto.arima()
命令,可以自动识别 ARIMA(p,d,q) 中各个阶数,非常方便。
library(fBasics)
basicStats(ret) # 对月度收益率的基本统计描述
t.test(ret) # 均值检验 H0: x=0
normalTest(ret,method='jb',na.rm=T) # 正态分布检验,JB-test
par(mfcol=c(1,1))
acf(ret, lag=12, plot=T) # 自相关系数作图
Box.test(ret,lag=6,type='Ljung') # m阶自相关系数为零的检验
library(forecast)
stock=auto.arima(ret,stationary = TRUE, seasonal = FALSE,
ic="aic")
# 拟合与预测
par(mfcol=c(2,1))
plot(stock$x, lty = 1, main = "上证综指收益率(%)")
lines(fitted(stock), lty = 2,lwd = 2, col = "red")
predict(stock, n.ahead=3)
plot(forecast(stock))