2016-11-24 28 views
0

BK帶通濾波器,我有時間Z系列在R和吉布斯現象

  Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec 
1922 -0.25108773 -0.27732553 -0.29703807 -0.30274000 -0.30323653 -0.28441682 -0.24106527 -0.18705071 -0.17440826 -0.17291725 -0.19116734 -0.21678948 
1923 -0.24487998 -0.26658925 -0.28613991 -0.29674346 -0.29335742 -0.28325761 -0.23326680 -0.18697904 -0.18443807 -0.18144226 -0.18190910 -0.21574376 
1924 -0.24465806 -0.27349425 -0.29925888 -0.30386766 -0.30250722 -0.27464960 -0.23390958 -0.19300616 -0.17910621 -0.17869576 -0.19611839 -0.20447324 
1925 -0.25326812 -0.27344637 -0.29352971 -0.30947682 -0.30872025 -0.27604449 -0.24065208 -0.19676031 -0.17172229 -0.18484153 -0.19542607 -0.21841577 
1926 -0.25214568 -0.27450911 -0.29438956 -0.30392114 -0.30619846 -0.29089168 -0.24829621 -0.20204202 -0.18621514 -0.18808172 -0.19708748 -0.22629595 
1927 -0.25107357 -0.27204514 -0.29494695 -0.30751442 -0.30800040 -0.28569694 -0.24655626 -0.19547608 -0.19018517 -0.18866641 -0.20132372 -0.22084811 
1928 -0.24733214 -0.27490388 -0.28780308 -0.30407576 -0.30857301 -0.28629658 -0.23872777 -0.19590465 -0.18437917 -0.18274289 -0.19936931 -0.22368973 
1929 -0.25531870 -0.27264628 -0.29418746 -0.30385231 -0.31022219 -0.27931003 -0.23404912 -0.19538227 -0.17226595 -0.18465123 -0.19072933 -0.22043396 
1930 -0.24735028 -0.27386782 -0.29193707 -0.29925459 -0.30039372 -0.28014958 -0.23551136 -0.19511701 -0.18006660 -0.18282789 -0.20113355 -0.22095253 
1931 -0.24903438 -0.27439043 -0.29219506 -0.30312159 -0.30557600 -0.28180333 -0.22676008 -0.19048014 -0.18982644 -0.18459638 -0.19550196 -0.22127202 
1932 -0.25870503 -0.27650825 -0.28521052 -0.30685609 -0.30896898 -0.28378619 -0.23614859 -0.18945699 -0.17575919 -0.17820312 -0.19620912 -0.21774873 
1933 -0.24187599 -0.25575287 -0.28325644 -0.29554461 -0.29018996 -0.27040369 -0.23514812 -0.19935749 -0.18732198 -0.18606057 -0.19327237 -0.22321366 
1934 -0.24793807 -0.26986056 -0.29217378 -0.30479126 -0.30199154 -0.27574924 -0.24097380 -0.18560708 -0.18643606 -0.18501770 -0.19375478 -0.22418002 
1935 -0.25587642 -0.27805131 -0.29239104 -0.30784907 -0.30459449 -0.28216514 -0.23839965 -0.20137460 -0.18619998 -0.18328896 -0.20121286 -0.22869388 
1936 -0.25322320 -0.28025116 -0.29713940 -0.30800346 -0.31177201 -0.28473251 -0.23552472 -0.20313945 -0.18251793 -0.18383941 -0.20554430 -0.23061875 
1937 -0.26268769 -0.28529769 -0.30230641 -0.31107806 -0.30183547 -0.28324508 -0.23840574 -0.19862786 -0.19297314 -0.19392849 -0.19603212 -0.22877177 
1938 -0.25445601 -0.28160871 -0.29837676 -0.29879519 -0.30328832 -0.28288226 -0.23577573 -0.19521124 -0.18393512 -0.19039895 -0.20537533 -0.21924241 
1939 -0.25180969 -0.28199995 -0.29601764 -0.30147945 -0.30372884 -0.27837795 -0.23720063 -0.19929773 -0.18770674 -0.19341142 -0.20753282 -0.22484697 
1940 -0.15145157 -0.16596690 -0.17572643 -0.18225920 -0.18823836 -0.17504012 -0.16019626 -0.12920340 -0.12369614 -0.12024704 -0.12891992 -0.14234080 
1941 -0.10045275 -0.11095497 -0.11585389 -0.11932455 -0.11976700 -0.11653216 -0.10259231 -0.08271703 -0.07621320 -0.07184160 -0.07284514 -0.07385666 
1942 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
1943 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
1944 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
1945 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
1946 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
1947 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
1948 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
1949 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
1950 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
1951 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
1952 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
1953 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
1954 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
1955 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
1956 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
1957 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
1958 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
1959 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
1960 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
1961 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 

enter image description here

,我想在9.7個月,16個月執行banpassfilter。當z是零的過濾器顯示仍有一些小循環我已經應用了bkfilter(包MFILTER)

enter image description here

然而,1942年後。在之前的研究中(Band-pass filter in R: weird behaviour at the end of time series),有人向我暗示,這種行爲可能是由吉布斯現象引起的。然後我已經改正了bkfunction這裏http://www.gla.ac.uk/media/media_219052_en.pdf

### Baxter-King filter 
modbkfilter <- function(x,pl=NULL,pu=NULL,nfix=NULL,typeBK=c("regular","modified"),type=c("fixed"),drift=FALSE) 
{ 
    if(is.null(drift)) drift <- FALSE 
    xname=deparse(substitute(x)) 
    type = match.arg(type) 

    if(is.null(type)) type <- "fixed" 

    if(is.ts(x)) 
     freq=frequency(x) 
    else 
     freq=1 

    if(is.null(pl)) 
    { 
     if(freq > 1) 
      pl=trunc(freq*1.5) 
     else 
      pl=2 
    } 

    if(is.null(pu)) 
     pu=trunc(freq*8) 

    b = 2*pi/pl 
    a = 2*pi/pu 

    n = length(x) 

    if(n<5) 
     warning("# of observations in Baxter-King filter < 5") 

    if(pu<=pl) 
     stop("pu must be larger than pl") 
    if(pl<2) 
    { 
     warning("in Baxter-King kfilter, pl less than 2 , reset to 2") 
     pl = 2 
    } 

    if(is.null(nfix)) 
     nfix = freq*3 

    if(nfix>=n/2) 
     stop("fixed lag length must be < n/2") 

    j = 1:(2*n) 
    if(typeBK=="regular") B = as.matrix(c((b-a)/pi,(sin(j*b)-sin(j*a))/(j*pi))) 
    if(typeBK=="modified") B = as.matrix(c(
             (b-a)/pi, 
             ((sin(j*b)-sin(j*a))/(j*pi)) * (sin((2*pi*j)/(2*nfix+1))/((2*pi*j)/(2*nfix+1))) 
             )) 

    AA = matrix(0,n,n) 

    if(type=="fixed") 
    { 
     bb = matrix(0,2*nfix+1,1) 
     bb[(nfix+1):(2*nfix+1)] = B[1:(nfix+1)] 
     bb[nfix:1] = B[2:(nfix+1)] 
     bb = bb-sum(bb)/(2*nfix+1) 

     for(i in (nfix+1):(n-nfix)) 
      AA[i,(i-nfix):(i+nfix)] = t(bb) 
    } 

    xo = x 
    x = as.matrix(x) 
    if(drift) 
     x = undrift(x) 

    x.cycle = AA%*%as.matrix(x) 
    x.cycle[c(1:nfix,(n-nfix+1):n)] = NA 
    x.trend = x-x.cycle 
    if(is.ts(xo)) 
    { 
     tsp.x = tsp(xo) 
     x.cycle=ts(x.cycle,star=tsp.x[1],frequency=tsp.x[3]) 
     x.trend=ts(x.trend,star=tsp.x[1],frequency=tsp.x[3]) 
     x=ts(x,star=tsp.x[1],frequency=tsp.x[3]) 
    } 
    res <- list(cycle=x.cycle,trend=x.trend,fmatrix=AA,title="Baxter-King Filter", 
       xname=xname,call=as.call(match.call()), 
       type=type,pl=pl,pu=pu,nfix=nfix,method="bkfilter",x=x) 

    return(structure(res,class="mFilter")) 
} 

描述然而,結果並沒有太大變化

enter image description here

任何幫助嗎?

回答

0

從方程三點紙:

equation 3

除數爲(2K + 1),但你已經使用(2*nfix + 1),由於j <= 1:(2*n),爲什麼不K <= 2*n(2*2*n + 1)在你的代碼呢?從我看到的nfix不等於2*n