R_misce

はじめに

• 講義や研究のために書いた小さなR関数を記そうと思います。(2021/09/15)
  • それにしてもこうして並べてみるとネーミングセンスのなさを思い知りますな。
  • 完全に無保証です。内容の変更も随時行い，それについて記さないことがあります。(2021/10/04)

• はじめに
• 関数とその使用例
• my.count.freqs()
• calcQp_plus()
• my.pcor(), my.pcor2()
• HPfilter()
• qHPfilter()
• qHPfilter()(Python版)

最終更新時間：2022年08月24日 08時15分20秒

関数とその使用例

 my.count.freqs()

• 階級下限以上，階級上限未満の度数を数え上げる関数です。最後の階級だけ，階級下限以上，階級上限以下の度数を数え上げます。(2021/09/15)

my.count.freqs=function(x,breaks){
  n=length(breaks)
  f=rep(0,times=n-1)
  for (i in 1:n-1){
     f[i]=length(subset(x,breaks[i]<=x & x<breaks[i+1]))
  }
  f[n-1]=f[n-1]+length(subset(x,x==breaks[n]))
  return(f)
}

> source("my.count.freqs.r")
> x=c(10,20,31,44,50,66,69,72,88,99)
> breaks=seq(0,100,by=20)
> print(my.count.freqs(x,breaks))
[1] 1 2 2 3 2

 calcQp_plus()

• 分位点を計算して返す関数です。quantile()関数の理解のために書きました。(2021/09/15)
  • ちなみに，quantile()関数のデフォルト・タイプは7です。

calcQp_plus = function(x,p,type=7){
x=sort(x)
n=length(x)
if (type==4) {h=n*p}
if (type==5) {h=n*p+1/2}
if (type==6) {h=(n+1)*p}
if (type==7) {h=(n-1)*p+1}
if (type==8) {h=(n+1/3)*p+1/3}
if (type==9) {h=(n+1/4)*p+3/8}
fl=floor(h)
ce=ceiling(h)
w=h-fl
Qp=(1-w)*x[fl]+w*x[ce]
return(Qp)
}

> x=rnorm(1000)
> source("calcQp_plus.r")
> calcQp_plus(x,0.2251,type=4); quantile(x,0.2251,type=4)
[1] -0.7112623
   22.51% 
-0.7112623 

 my.pcor(), my.pcor2()

• 偏相関係数計算のための関数です。(2021/09/15)

my.pcor=function(x,y,z){
iota=matrix(1,nrow=nrow(z)); Z=cbind(iota,z)
res_x=x-Z%*%solve(crossprod(Z),t(Z)%*%x)
res_y=y-Z%*%solve(crossprod(Z),t(Z)%*%y)
pc=cor(res_x,res_y)
return(pc)
}

my.pcor2=function(x,y,z){
m.x=lm(x~z)
m.y=lm(y~z)
pc=cor(m.x$residuals,m.y$residuals)
return(pc)
}

 HPfilter()

• ホドリック・プレスコット・フィルターの関数です。(2021/09/15)

HPfilter <- function(y,lambda){
  n <- length(y)
  I <- diag(rep(1,times=n))
  D<-diff(I,differences=2)
  xhat <- solve(I+lambda*crossprod(D),y)
  return(xhat)
}

• HPfilter.r(53)
• HPfilter_r_smpl.r(55)
• dataSS.txt(57)

 qHPfilter()

• HPフィルターの分位点版です。論文，投稿中です。上手くいくことを祈っています。(2021/09/15)←Empirical Economics誌に掲載されました(リンク)。(2022/08/24)
  • 大雑把に言って，残差のうち100tau%が負の値になるようにスムーズなトレンド線が引かれます。(下の累積相対度数のグラフ参照。)(2021/09/17)
  • CVXRパッケージ(リンク)が必要です。

library(CVXR)

qHPfilter=function(y,lambda,tau){
# y: Tx1 matrix
# lambda: positive real number
# tau: real number belongs to (0,1)
T=nrow(y)
D=diff(diag(T),differences=2)
x=Variable(T)
p=Problem( Minimize( 0.5*p_norm(y-x,1)+(tau-0.5)*sum(y-x)+lambda*sum_squares(D%*%x) ) )
result=solve(p)
xhat=result$getValue(x)
return(xhat)
}

• qHPfilter.r(57)
• qHPfilter_smpl.r(59)
• dataSS.txt(57)

 qHPfilter()(Python版)

• Python版を試作してみました。(2021/10/04)
  • (NumPyのほか)CVXPY(リンク)が必要です。(2021/10/04)
  • Pythonは勉強し始めたばかりで今一つよく分かりません(^^;。(2021/10/04)

import numpy as np
import cvxpy as cp

def qHPfilter(y,lam,tau):
# y: Tx1 matrix
# lam: positive real number
# tau: real number belongs to (0,1)
    T=y.size
    D=np.diff(np.eye(T),n=2).T
    x=cp.Variable(T)
    p=0.5*cp.norm(y-x,1)+(tau-0.5)*cp.sum(y-x)+lam*cp.sum_squares(D@x)
    prob=cp.Problem(cp.Minimize(p))
    prob.solve(verbose=True)
    return x.value

year=np.loadtxt('./dataSS1.txt')
y=np.loadtxt('./dataSS2.txt')

tau09=0.9
tau05=0.5
tau01=0.1
lam=10**4
xhat09=qHPfilter(y,lam,tau09)
xhat05=qHPfilter(y,lam,tau05)
xhat01=qHPfilter(y,lam,tau01)

import matplotlib.pyplot as plt

plt.figure(figsize=(6, 6))
plt.plot(year,y,'k-',linewidth=1.0)
plt.plot(year,np.array(xhat09),'r-',linewidth=1.0)
plt.plot(year,np.array(xhat05),'g-',linewidth=1.0)
plt.plot(year,np.array(xhat01),'b-',linewidth=1.0)
plt.show()
#plt.savefig("SST_py.png")
#plt.savefig("SST_py.pdf")

• qHPfilter_smpl.py(45)(2021/10/04)
• dataSS1.txt(40)
• dataSS2.txt(44)