R Time Series Analysis 時系列解析(2) (decompose)

#decompose関数を分解して理解する
#decompose関数は系列データを季節成分、トレンド、残りに分解する関数
#データは経済産業省より取得した「総合原指数【月次】出荷（平成17年＝100.0）」
edit(decompose)

dat <- read.csv("http://www.meti.go.jp/statistics/tyo/iip/result/h2afdldj/csv/ha2zom3j.csv",
                skip=2,header=TRUE) #skip=2は上位2行を省く意味
dat_ts <- ts(dat[,4],start=c(2003,1),frequency=12)

l <- length(dat_ts)
f <- frequency(dat_ts)
type <- "additive"
#季節成分(sesonal),トレンド(trend),残り(random)をどのように分解するか。
#"additive"なら和型、"multiplicative"なら積型
filt <- NULL

if (f <= 1 || length(na.omit(dat_ts)) < 2 * f) {
  stop("time series has no or less than 2 periods")  
}
# ||は「または」の意味

#filter(移動平均時の重み)の設定
if (is.null(filt)) {
  #filterが設定されていない場合は下記の通り設定設定
    filt <- if (!f%%2) {#f%%2はfを2で割った余り
      #fが2で割り切れない場合
      c(0.5, rep(1, f - 1), 0.5)/f #1周期分+1で中央移動平均をとる。両端は0.5/f,中央は1/fとする。
    } else {
      #fが2で割り切れる場合
      rep(1, f)/f #中央移動平均の重みは一律1/f
    }
}

#trendの抽出
    #設定した重み(filter)でdat_tsの中心移動平均をとり、"trend"とする
    trend <- filter(dat_ts, filt)

#seasonの抽出
    season <- if (type == "additive") {
      dat_ts - trend
    } else {
      dat_ts/trend
    }

periods <- l%/%f　#周期数を算出
index <- seq.int(1L, l, by = f) - 1L #12カ月毎の整数列を生成
figure <- numeric(f) #12個の数値の入れ物

#月ごとに
for (i in 1L:f) { 
  figure[i] <- mean(season[index + i], na.rm = TRUE) #各月毎の平均をNAを除いて計算
}

figure <- if (type == "additive") {
  figure - mean(figure)  #"addtive"の場合
} else {
  figure/mean(figure)    #"multiplicative"の場合
}

seasonal <- ts(rep(figure, periods + 1)[seq_len(l)],　#各付毎の平均をseasonalとする
               start = start(dat_ts), frequency = f)

#原系列(dat_ts),季節成分(sesonal),トレンド(trend),残り(random)を構造化
structure(list(x = dat_ts, seasonal = seasonal, trend = trend, 
               random = if (type == "additive") dat_ts - seasonal - trend else dat_ts/seasonal/trend, 
               figure = figure, type = type), class = "decomposed.ts")