!!!はじめに

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

{{outline}}

{{lastmodified}}

!!!関数とその使用例

!!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)
 }

*{{ref HPfilter.r}}
*{{ref HPfilter_r_smpl.r}}
*{{ref dataSS.txt}}

{{ref_image HPfilter.png,w300,h300}}

!!qHPfilter()

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

 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)
 }

*{{ref qHPfilter.r}}
*{{ref qHPfilter_smpl.r}}
*{{ref dataSS.txt}}

{{ref_image qHPfilter.png,w300,h300}}
{{ref_image qHP2.png,w330,h330}}

!!qHPfilter()(Python版)

*Python版を試作してみました。(2021/10/04)
**(NumPyのほか)CVXPY([リンク|https://anaconda.org/conda-forge/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")

*{{ref qHPfilter_smpl.py}}(2021/10/04)
*{{ref dataSS1.txt}}
*{{ref dataSS2.txt}}

{{ref_image SST_py.png,w300,h300}}
----

{{attach}}