R prg of Regression Kink with an Unknown Threshold
对R的版本要求不高,只需要在百度R就下载base版就可以跑了,详细的data和prg可以直接戳 download
由于了解的门限节点模型,只需要其中一个growth.r的程序就可以完成了。R语言的基本操作可以百度or找一本基础的书比如《R语言实战》很基础,也很实用。
源程序
growth.r
This is an R file
It creates the empirical work reported in Bruce Hansen “Regression Kink with an Unknown Threshold”
It uses the data file “usdata.txt”, extracted from the Reinhart-Rogoff file “41_data.xlsx”
For other applications, load in your own data into the vectors y, x, and matrix z
Include an intercept in z
The vector gammas is the set of thresholds to search
The vector dx is for displaying the regression estiamtes
—Load in Data
growth <- read.table(“c:/Users/Tinho/Desktop/文献/2015-Regression Kink with an Unknown Threshold/2015-Regression Kink with an Unknown Threshold/cthresh/usdata.txt”,header=TRUE)
n = nrow(growth)
year = growth[2:n,1]
gdp = growth[2:n,3]
gdp1 = growth[1:(n-1),3]
debt1= growth[1:(n-1),2]
—Time-Series Data Plots, Figure 1ab in paper
windows()
plot(year,gdp,type=”l”,ylab=”GDP Growth Rate”)
savePlot(file=”fig1a.eps”,type=”eps”,dev.cur())
windows()
plot(year,debt1,type=”l”,ylab=”Debt/GDP”)
savePlot(file=”fig1b.eps”,type=”eps”,dev.cur())
— Define variables
y = gdp
x = debt1
n = length(y)
z = cbind(gdp1,matrix(1,n,1))
—Controls
gammas = seq(10,70,by=0.1) # Grid on Threshold parameter for estimation
dx = seq(0,121,by=1) # Grid on regression function for display
level = 0.90 # For confidence sets
boot = 10000 # Number of bootstrap replications
Ceps = c(0.5,1,2,4) # For numerical delta method bootstrap
############################################
— Some useful functions
reg <- function(X,y) {
X <- qr(X)
as.matrix(qr.coef(X,y))
}
pos.part <- function(x) x(x>0)
neg.part <- function(x) x(x<0)
pt1 = proc.time()
— Linear Model
x0 = cbind(debt1,z)
k0 = ncol(x0)
kz = ncol(z)
x00 = solve(crossprod(x0))
bols = x00%%crossprod(x0,y)
e0 = y - x0%%bols
sse0 = sum(e0^2)
sigols = sse0/n
v0 = x00%%crossprod(x0matrix(e0,n,k0))%%x00(n/(n-k0))
seols = as.matrix(sqrt(diag(v0)))
— Threshold Model
grid = length(gammas)
rd = length(dx)
sse = matrix(0,grid,1)
k = kz + 3
for (j in 1:grid) {
gamj=gammas[j]
x1 = cbind(neg.part(x-gamj),pos.part(x-gamj),z)
e1 = y - x1%%reg(x1,y)
sse[j] = sum(e1^2)
}
gi = which.min(sse)
gammahat = gammas[gi]
ssemin = sse[gi]
x1 = cbind(neg.part(x-gammahat),pos.part(x-gammahat),z)
bt = reg(x1,y)
et = y - x1%% bt
hg = - (x
x2 = cbind(x1,hg)
hg2 = crossprod(cbind((x
xx2 = matrix(0,k,k)
xx2[1:2,k]=hg2
xx2[k,1:2]=t(hg2)
xxi = solve(crossprod(x2) + xx2)
v = xxi%
betahat = rbind(bt,gammahat)
se = as.matrix(sqrt(diag(v)))
sig = sum(et^2)/n
wt = n(sse0-ssemin)/ssemin
wg = n*(sse-ssemin)/ssemin
— Plot Least Squares Criterion, Figure 3 in paper
windows()
plot(gammas,sse,type=”l”,ylab=”Least Squares Criterion”,xlab=”Threshold Parameter”)
savePlot(file=”fig3.eps”,type=”eps”,dev.cur())
—Regression Estimate
G = cbind(neg.part(dx-gammahat),pos.part(dx-gammahat),matrix(1,rd,1)%%colMeans(z))
yf = G%%bt
— Bootstrap & Testing
waldb = matrix(0,grid,boot)
sseb = matrix(0,grid,boot)
betab = array(0,c(grid,k-1,boot))
u = matrix(rnorm(nboot),n,boot)
eb = matrix(e0,n,boot)u
yb = matrix(x1%%bt,n,boot) + matrix(et,n,boot)u
eb0 = eb - x0%%reg(x0,eb)
bsse0 = colSums(eb0^2)
for (j in 1:grid) {
gamj = gammas[j]
x2 = cbind(neg.part(x-gamj),pos.part(x-gamj),z)
eb0 = eb - x2%%reg(x2,eb)
bsse = colSums(eb0^2)
waldb[j,] = n(bsse0-bsse)/bsse
bb = reg(x2,yb)
eb1 = yb - x2%%bb
sseb[j,] = colSums(eb1^2)
betab[j,,] = bb
}
— Multiplier Bootstrap test for Threshold
wb = apply(waldb,2,max)
pv = mean(wb > matrix(wt,boot,1))
crit = quantile(wb,probs=level)
— Threshold Regression Estimates
gib = apply(sseb,2,which.min)
gamb = gammas[gib]
— Symmetric Percentile Confidence Interval Construction
betahatb = matrix(0,k-1,boot)
ci = matrix(0,k,2)
for (j in 1:(k-1)){
bj = diag(betab[gib,j,])-bt[j]
qj = quantile(abs(bj),probs=level)
ci[j,1] = bt[j] - qj
ci[j,2] = bt[j] + qj
betahatb[j,] = t(bj)
}
— Confidence Interval for Threshold
sseb0 = colSums((yb - x1%%reg(x1,yb))^2)
sseminb = apply(sseb,2,min)
wgb = n(sseb0-sseminb)/sseminb
qa = qchisq(level,1)
wia = (wg > qa)
cga = c(gammas[which.min(wia)],gammas[grid+1-which.min(rev(wia))])
qb = quantile(wgb,probs=level)
wib = (wg > qb)
cgb = c(gammas[which.min(wib)],gammas[grid+1-which.min(rev(wib))])
ci[k,] = cgb
— Threshold Regression Confidence Intervals
mdx = t(matrix(dx,rd,boot))-gammahat
thetab = t(G%%betahatb)
cn = length(Ceps)
yf1 = matrix(0,rd,cn)
yf2 = matrix(0,rd,cn)
for (j in 1:cn) {
eps = Ceps[j]
h = matrix(gamb-gammahat,boot,rd)eps
thetaq = thetab + ((neg.part(mdx-h)-neg.part(mdx))bt[1] + (pos.part(mdx-h)-pos.part(mdx))bt[2])/eps
qf = apply(abs(thetaq),2,quantile,probs=level)
yf1[,j] = yf - qf
yf2[,j] = yf + qf
}
pt2 = proc.time()
sink(“growth.out”)
cat(“Linear Model, coefficients and error variance”,”\n”)
print(cbind(bols,seols),digits=2)
cat(“Error variance”,”\n”)
print(sigols,digits=4)
cat(“Wald Test for Threshold, p-value, & critical value”,”\n”)
print(c(wt,pv,crit,level),digits=3)
cat(“Threshold Model Estimates, s.e.’s, and Bootstrap confidence intervals”,”\n”)
print(cbind(betahat,se,ci),digits=2)
cat(“Error variance”,”\n”)
print(sig,digits=4)
cat(“Bootstrap Critical value for threshold parameter interval”,qb,”\n”)
cat(“Computation Time: replications and seconds”,”\n”)
print(c(boot,pt2[3]-pt1[3]),digits=3)
sink()
— Plot Confidence Interval Construction for Threshold (Figure 4)
windows()
plot(gammas,wg,type=”l”,ylab=”Threshold F Statistic”,xlab=”Threshold Parameter”)
cr1 = matrix(qa,grid,1)
cr2 = matrix(qb,grid,1)
lines(gammas,cr1,lty=”dashed”,col=”blue”)
lines(gammas,cr2,lty=”dotted”,col=”red”)
arrows(cga[1],qa,cga[1],0,lty=”dashed”,col=”blue”)
arrows(cga[2],qa,cga[2],0,lty=”dashed”,col=”blue”)
arrows(cgb[1],qb,cgb[1],0,lty=”dotted”,col=”red”)
arrows(cgb[2],qb,cgb[2],0,lty=”dotted”,col=”red”)
legend(“topright”,c(“Bootstrap Critical Value”,”Asymptotic Critical Value”),lty=c(“dotted”,”dashed”),col=c(“red”,”blue”))
savePlot(file=”fig4.eps”,type=”eps”,dev.cur())
—Scatter plot, regression line, and confidence intervals (Figure 2)
windows()
plot(debt1,gdp,ylab=”GDP Growth Rate”,xlab=”Debt/GDP”)
lines(dx,yf)
lines(dx,yf1[,2],lty=”dashed”,col=”blue”)
lines(dx,yf2[,2],lty=”dashed”,col=”blue”)
yk = colMeans(z)%*%bt[3:(2+kz)]
points(gammahat,yk,col=”red”,bg=”red”,pch=22)
savePlot(file=”fig2.eps”,type=”eps”,dev.cur())
windows()
plot(dx,yf,ylab=”GDP Growth Rate”,xlab=”Debt/GDP”,type=”l”,ylim=c(-8,8),yaxt=”n”)
axis(2,at=c(-8,-4,0,4,8))
lines(dx,yf1[,1],lty=”dashed”,col=”orange”)
lines(dx,yf2[,1],lty=”dashed”,col=”orange”)
lines(dx,yf1[,2],lty=”dotted”,col=”red”)
lines(dx,yf2[,2],lty=”dotted”,col=”red”)
lines(dx,yf1[,3],lty=”dashed”,col=”blue”)
lines(dx,yf2[,3],lty=”dashed”,col=”blue”)
lines(dx,yf1[,4],lty=”dotted”,col=”black”)
lines(dx,yf2[,4],lty=”dotted”,col=”black”)
points(gammahat,yk,col=”red”,bg=”red”,pch=22)
leg1 = bquote(paste(epsilon[n]==.(Ceps[1]),n^-.5))
leg2 = bquote(paste(epsilon[n]==.(Ceps[2]),n^-.5))
leg3 = bquote(paste(epsilon[n]==.(Ceps[3]),n^-.5))
leg4 = bquote(paste(epsilon[n]==.(Ceps[4]),n^-.5))
leg_text = c(as.expression(leg1),as.expression(leg2),as.expression(leg3),as.expression(leg4))
legend(“bottomleft”,leg_text,lty=c(“dashed”,”dotted”,”dashed”,”dotted”),col=c(“orange”,”red”,”blue”,”black”))
savePlot(file=”fig7.eps”,type=”eps”,dev.cur())