x12simldoc92
2023 年 1 月 29 日

1

X12SIML を使った季節調整の実際

佐藤整尚 (東京大学)

1.1

X12SIML (Ver.92)• SIML 法を使った季節調整プログラムにセンサス局の X12ARIMA のような異常値処理を加えたもの• R 上で動作する。
• AO,LS,RAMP などの自動検出機能あり。
• ユーザー定義の回帰変数を入れることも可能• source(“x12siml92.R”) でインストールする。その他ライブラリ等は必要ない。
[1]: source("x12siml92.R")

1.2

用法“x12siml”<-

function(data, reg = NULL, trend = 2, ilog = 0, frequency = 4, start= c(1994, 1), iplot = T, sorder = 1, mtype=1,

pb=2,pa=2, …)• reg のところにはユーザー定義の回帰変数を入れることができる。
• trend, sorder には平滑化の度合いを入れる。大きいほうがより柔軟になる。目安として trendにはデータ数の１０％ぐらい、sorder は１～３を指定。

• pb で指定された回数だけ、データの始めの１周期を繰り返す。pa で指定された回数だけ、
データ最後の１周期を繰り返す。

• mtype で SIML フィルターのタイプを選ぶ。=1 (前向き) or =2 (後向き).
• … のところに以下の outlier を示す変数を入れることができる。
• 構造変化・異常値処理– AO（加法的外れ値）– LS（レベルシフト）– TC（一時的変化）1

– RAMP（傾斜変化）– VAT– Holiday– Leap Year
– TCRAMP(RAMP+TC)– TCRAMP1(RAMP+(１時点空く)+TC)• outlier 変数• AO > x12siml(data,start=c(1993,1),frequency=4,ao=c(2001,2))• 2001

年 第 ２ 四 半 期 にAO

を セ ッ ト>x12siml(data,start=c(1993,1),fre-quency=4,aot=c(2001,2))• 2001 年第２四半期に AO をセットし、トレンドに入れる。
• LS > x12siml(data,start=c(1993,1),frequency=4,ls=c(2001,2))• 2001 年第２四半期に LS をセット
• RAMP > x12siml(data,start=c(1993,1),frequency=4,rp=c(2001,2,2001,4))• 2001 年第２四半期から第４四半期にかけて RAMP をセット
• TLS > x12siml(data,start=c(1993,1),frequency=4,tls=c(2001,2,2001,4))• 2001 年第２四半期から第４四半期にかけて一時的レベルシフト (TLS) をセット
• TC >

x12siml(data,start=c(1993,1),frequency=4,tc=c(2001,2))• 2001 年第２四半期に TC をセット
• VAT >

x12siml(data,start=c(1993,1),frequency=4,vat=c(2001,1))• 消 費 税 ダ ミ ー （指 定 し た 時 点 が 1、 次 の 時 点 が-1 と い う ダ ミ ー)>

x12siml(data,start=c(1993,1),frequency=4,vatt=c(2001,1))• 消費税ダミー（指定した時点が 1、次の時点が-1 というダミー)、ただし、トレンドに入れる。
• Holiday >x12siml(data,start=c(1993,1),frequency=4,hol=c(2001,3,2003,4))• AO と同じダミーで複数時点指定可能。推定結果は AO と違って、季節性に含める。
• LeapYear >x12siml(data,start=c(1993,1),frequency=4, ly=c(1996,1))• 指定された時点とそれ以降 4*frequency• SLS >

ごとに AO をセット。

x12siml(data,start=c(1993,1),frequency=4, sls=c(2003,4))• 指定された時点に季節性レベルシフトをセット。
• TCRAMP > x12siml(data,start=c(1993,1),frequency=4,tcrp=c(2001,2,2001,4,2003,3))• 2001 年第２四半期から 2002 年第３四半期にかけて TCRAMP をセット。2001 年第２四半期
から 2001 年第 4 四半期は RAMP、2001 年第 4 四半期から 2003 年第３四半期まで TC• TCRAMP1 > x12siml(data,start=c(1993,1),frequency=4,tcrp1=c(2001,2,2001,4,2003,3))• 2001 年第２四半期から 2002 年第３四半期にかけて TCRAMP1 をセット。2001 年第２四半
期から 2001 年第 4 四半期は RAMP、2002 年第 1 四半期から 2003 年第３四半期まで TC

[ ]:

2

[2]: TCRATE <- 0.8^2start <- c(1993,1)
frequency <- 4
reg2 <- list(ao=c(2001,2),ls=c(2001,2),rp=c(2001,2,2001,4),tc=c(2001,2), vat=c(2001,1),hol=c(2001,3,2003,4), sls=c(2003,4),
ly=c(1996,1),

tcrp=c(2001,2,2001,4,2003,3),tls=c(2001,2,2001,4),

tcrp1=c(2001,2,2001,4,2003,3))[5]: data <- rnorm(100)n <- length(data)nn <- length(reg2)
dimnames.reg <- NULL
reg <- NULL
if(nn > 0) {for(i in seq(nn)) {switch(names(reg2[i]),"ls"={z <- rep(0, n)
tt <- sum((reg2[[i]] - start) * c(frequency,␣
,→

1)
) + 1
z[1:(tt - 1)] <- -1
z <- z+1
reg <- cbind(reg, z)},"tc"= {
tcrate <- TCRATE #0.69999999999999996^(12/,→

frequency)z <- rep(0, n)
tt <- sum((reg2[[i]] - start) * c(frequency,␣

,→

1)
) + 1
z[tt:n] <- tcrate^(0:(n - tt))
reg <- cbind(reg, z)},

3

"ao"= {z <- rep(0, n)
tt <- sum((reg2[[i]] - start) * c(frequency,␣
,→

1)
) + 1
z[tt] <- 1
reg <- cbind(reg, z)},"hol"= {z <- rep(0, n)
nnn2 <- length(reg2[[i]])/2
for(j in seq(nnn2)) {tt <- sum((reg2[[i]][j*2-(1:,→0)]-start)*c(frequency,1)) + 1
z[tt] <- 1}
reg <- cbind(reg, z)},"sls"= {z <- rep(0, n)
nyy <- floor(n/frequency)+1
yyy <- c(rep(start[1],frequency-start[2]+1),t(matrix(rep((start[1]+1):(start[1]+nyy),frequency),
nrow=nyy)))[1:n]mmm <- c((start[2]):frequency,rep(1:frequency,nyy))[1:n]sty <- reg2[[i]][1]stm <- reg2[[i]][2]z[yyy >= sty] <- -1/(frequency-1)z[yyy >= sty & mmm==stm] <- 1
reg <- cbind(reg, z)},"ly"= {z <- rep(0, n)if(length(reg2[[i]])==1)
tt <- reg2[[i]]+(0:n)*(frequency*4)

4

else
tt <- sum((reg2[[i]] - start) * c(frequency,␣
,→1)) + 1 + (0:n)*(frequency*4)z[tt] <- 1
z <- z[1:n]reg <- cbind(reg, z)},"vat"={z <- rep(0, n)
tt <- sum((reg2[[i]] - start) * c(frequency,␣

,→

1)
) + 1
z[tt] <- 1z[tt + 1] <- -1
reg <- cbind(reg, z)},"rp"= {z <- rep(0, n)
tt1 <- sum((reg2[[i]][1:2] - start) * c(frequency, 1)) + 1
tt2 <- sum((reg2[[i]][3:4] - start) * c(frequency, 1)) + 1
z[1:(tt1)] <- -1z[(tt1 + 1):(tt2 - 1)] <- (((tt1 + 1):(tt2 -␣,→

1)) - tt1)/(tt2 - tt1) - 1
z <- z+1
reg <- cbind(reg, z)},"tls"= {z <- rep(0, n)
tt1 <- sum((reg2[[i]][1:2] - start) * c(frequency, 1)) + 1
tt2 <- sum((reg2[[i]][3:4] - start) * c(frequency, 1)) + 1

5z[(tt1):tt2] <- 1
reg <- cbind(reg, z)},"tcrp"= {z <- rep(0, n)
tt1 <- sum((reg2[[i]][1:2] - start) * c(frequency, 1)) + 1
tt2 <- sum((reg2[[i]][3:4] - start) * c(frequency, 1)) + 1
tt3 <- sum((reg2[[i]][5:6] - start) * c(frequency, 1)) + 1
tcrate <- TCRATE#z[1:(tt1)] <- -1z[(tt1 + 1):(tt2 - 1)] <- (((tt1 + 1):(tt2 -␣,→

1)) - tt1)/(tt2 - tt1)z[(tt2):tt3] <- tcrate^(0:(tt3 - tt2))z[tt3:n] <- z[tt3]reg <- cbind(reg, z)},"tcrp1"= {z <- rep(0, n)
tt1 <- sum((reg2[[i]][1:2] - start) * c(frequency, 1)) + 1
tt2 <- sum((reg2[[i]][3:4] - start) * c(frequency, 1)) + 1
tt3 <- sum((reg2[[i]][5:6] - start) * c(frequency, 1)) + 1
tcrate <- TCRATE#0.80000000000000004##^(12/frequency)z[1:(tt1)] <- -1z[(tt1 + 1):(tt2 - 1)] <- (((tt1 + 1):(tt2 -␣,→

1)) - tt1)/(tt2 - tt1)
z[tt2] <- 1z[(tt2 + 1):tt3] <- tcrate^(0:(tt3 - tt2 -␣

,→

1))

6

z[tt3:n] <- z[tt3]reg <- cbind(reg, z)},{
nn <- nn-1cat(paste("Warning:",names(reg2[i]), "is not␣,→

supported\n"))reg2[i] <- NULL})}}
dimnames(reg) <- list(NULL,names(reg2))[6]: par(mfrow=c(3,1))
plot.ts(reg[,"ao"],main="AO")
plot.ts(reg[,"ls"],main="LS")plot.ts(reg[,"rp"],main="RAMP")7[11]: par(mfrow=c(3,1))
plot.ts(reg[,"tc"],main="TC")plot.ts(reg[,"vat"],main="VAT")
plot.ts(reg[,"hol"],main="Holiday")8[18]: par(mfrow=c(3,1))
plot.ts(reg[,"ly"],main="LeapYear")plot.ts(reg[,"tcrp"],main="TCRAMP")
plot.ts(reg[,"tcrp1"],main="TCRAMP1")9[6]: par(mfrow=c(3,1))
plot.ts(reg[,"tls"],main="TLS")
plot.ts(reg[,"sls"],main="SLS")

10

1.3

季節周波数
• 季節周期を s とする。（つまり、月次：s = 12 、四半期：s = 4 ）• 季節性が強いと思われる季節周波数 fs は
1 2
[s/2]
fs = ( , , · · · ,
)
s s
s
1 2
6
f12 = ( , , · · · , )
12 12
12
1 2 3
f7 = ( , , )
7 7 7
1 2
f4 = ( , )4 4• 揺らぎの幅（sorder）を考慮した場合の季節インデックスを
inds = (2 ∗ fs ∗ n − sorder) : (2 ∗ fs ∗ n + sorder)(≤ n)11

ind4 = (n/2 − sorder) : (n/2 + sorder), (n − sorder) : nとする。ただし、Ｘ軸の１～ｎが周波数 0~1/2 に対応していると仮定している。

[2]: h <- 0
n1 <- 100
frequency <- 8sorder <- 2
period <- frequency
if(F) {
sid <- ceiling(2*(n1+h)/period)+(-sorder):(sorder)
if(period==4) { sid <- c(sid,n1-(( max(sorder+1-h,0) ):0))}if(period==12) {sid <- c(sid,ceiling(4*(n1+h)/period)+(-sorder):(sorder))sid <- c(sid,ceiling(6*(n1+h)/period)+(-sorder):(sorder))sid <- c(sid,ceiling(8*(n1+h)/period)+(-sorder):(sorder))sid <- c(sid,ceiling(10*(n1+h)/period)+(-sorder):(sorder))sid <- c(sid,n1-(( max(sorder+1-h,0) ):0))}if(period==7) {sid <- c(sid,ceiling(4*(n1+h)/period)+(-sorder):(sorder))sid <- c(sid,ceiling(6*(n1+h)/period)+(-sorder):(sorder))## sid <- c(sid,n1-(( max(sorder+1-h,0) ):0))}}sid <- ceiling(2*(n1+h)/period)+(-sorder):(sorder)
if(period > 4) {for(i in 2:(ceiling(period/2)-1)) {sid <- c(sid,ceiling(2*i*(n1+h)/period)+(-sorder):(sorder))}}if(period %% 2 == 0)

sid <- c(sid,n1-(( max(sorder+1-h,0) ):0))

sfreq <- rep(0,n1)
sfreq[sid] <- 1
plot.ts(sfreq,type="h")abline(h=0)

12

1.4

回帰 (regsiml)• A) データを Kn 変換したものを Z、回帰変数も Kn 変換を行い、Z (d) を得る。
z = P C −1 y, Zd = P C −1 W• B) Ｚ系列のうち、季節インデックスに対応する部分のみを抜きだした Zs を作る。
z s = Js z
Js は季節周波数の周りだけを取り出す変換• C) Ｚのうちに季節インデックスに対応する部分以外を抜き出したものを Z−s とする。
z−s = J−s z
J−s は季節周波数の周り以外を取り出す変換13• D) Z−s のうち始めのｍ個（m は引数 trend で指定) を抜きだした Zt を作る。
zt = Jm z
Jm は始めのｍ個だけを取り出す変換• E) Z (d) に対しても同様にして、Zs ,Z−s ,Zt
(d)

(d)

(d)

を作る。

Zds = Js Zd, Zd−s = J−s Zd, Zdt = Jm Zd• F) 回帰Ａ：Zt を、Zt(d)に回帰する (切片 0)。

Ct = (Zd0t Zdt )−1 Zd0t zt• G) 回帰Ｂ：Z−s を、Z−s に回帰する (切片 0)。
(d)C−s = (Zd0−s Zd−s )−1 Zd0−s z−s• H) ダミー変数の内、AO,VAT,HOL,LY のいずれかを含む場合は、Ｆ値を比較して、回帰Ａ，
Ｂのうち良いほうを選ぶ。それ以外は回帰 A を行う。

•

I) その回帰での回帰係数 coef と残差

res を得る。

0
C = Ct orC−s , res = Jm0 zt − Zdt Ct orJ−s
z−s − Zd−s C−s• J) Zs(d)に coef をかけて Zs から引いたものを Zs∗ とする。

zs ∗ = zs − Zds C• K) res のうち始めのｍ個を Zt∗ とする。（残りの n1-m 個は 0 とする。）
zt ∗ = Jm res• L) Zs∗ と Zt∗ を逆変換したものを、季節性Ｓ、トレンド Tr とする。
0
S = CP 0 Js0 zs ∗, T r = CP 0 Jm
zt ∗• M) 回帰 A の AIC を AIC.a、回帰 B の AIC を AIC.b とする。( AIC を求めるとき、分散の Log に n1 を掛けるーーーｍではなく）[8]: ## regsiml
mat <- siml.mat(n,m1,type=mtype)## type=1 (前向き), 2 (後ろ向き)
n1 <- dim(mat$K)[1]14

## A)
z.y <- mat$K %*% c(data)z.d <- mat$K %*% reg
z.s <- z.y## B)z.s[-sid] <- 0## C)
z.y[sid] <- 0## D)E)F)
zz <- lsfit(z.d[1:m1,],z.y[1:m1],inter=F)
vvv2 <- mean((zz$res)^2)##G)
zz.tmp <- lsfit(z.d[-sid,],z.y[-sid],inter=F)
vvv3 <- mean((zz.tmp$res)^2)
if( "F of zz" < "F of zz.tmp" and (ao or vat or ly or hol) ) zz <- zz.tmp## K)
res <- c(zz$res,rep(0,n1-m1))## J)
z.s[sid] <- z.s[sid]-c(z.d %*% zz$coef)[sid]## L)trend <- mat$inv %*% resseasonal <- mat$inv %*% z.s## M)para <- log(c(vvv2,vvv3))*n1+2*(k+1)
Error in makeMat21(n, m, type = type): object 'm1' not foundTraceback:1. siml.mat(n, m1, type = mtype)2. makeMat21(n, m, type = type)

15

1.5

結果

(x12siml)-

x12siml では regsiml の結果を加工して、図示をする。+ ダミー変数を含む回帰部分の内、“ao”,“vat”,“hol”,“ly” 以外の部分は、トレンドに含まれるようにする。+に含まれるようにする。+“hol”,“ly”

は季節性“ao”,“vat” はノイズに含まれるようにする。

[ ]: zz <- regsiml(data, reg, m1 = trend, log = ilog, period =frequency, sorder = sorder,
mtype=mtype, pb=pb, pa=pa)
if(ilog > 0)data <- log(data)z.trend <- zz$trend
if(!is.null(reg)) {zz.dumm <- seq(ncol(zz$trade))zz.dumm <- zz.dumm[dimnames(zz$trade)[[2]] != "ao" &dimnames(zz$trade)[[2]] != "vat" &dimnames(zz$trade)[[2]] != "hol" &dimnames(zz$trade)[[2]] != "ly"]if(length(zz.dumm) >0) {z.trend <- z.trend + c(apply(zz$trade[,zz.
,→dumm,drop=F],1,sum)) }}zz.seasonal <- zz$seasonal
if(any(dimnames(reg)[[2]] == "ly"))zz.seasonal <- zz.seasonal+zz$trade[,"ly"]if(any(dimnames(reg)[[2]] == "hol")) {ididi <- seq(nn0)[dimnames(reg)[[2]] == "hol"]for(ijij in ididi) {zz.seasonal <- zz.seasonal+zz$trade[,ijij]}}
z.adj <- data - zz.seasonalzz.noise <- data-(zz.seasonal+z.trend)

16

1.6

異常値自動探索outlier(data, start=c(1998,1), frequency=12,type=“ao”,tt=c(1999,1,2009,12))• AO を 1999 年 1 月～2009 年 12 月の間で自動探索。AIC がある程度以上（デフォルト４）よ
くなる点のうち、もっとも、改善する点 1 点のみを探索。

• 複 数 検 出 す る 場 合 は、>

outlier(data,

start=c(1998,1),

fre-

quency=12,type=“ao”,tt=c(2001,4,2009,12), ao=c(2001,3))–（期間がだぶらないように注意する。）• type には“ao”, “ls”, “tc”, “rp”

を指定可

• TKRAMP、TCRAMP の自動探索
tkramp(data, start=c(1996,1), frequency=4,tt=c(2008,1,2010,12))• tkramp を 2008 年 1 月～2010 年 12 月の間で自動探索。（長い期間を指定すると、時間がかかる）• TKRAMP とはダブル RAMP（高岡 RAMP、RAMP+RAMP）• AIC がある程度以上（デフォルト４）よくなる期間のうち、もっとも、改善する期間 1 期間のみを探索。

• ほ か の 異 常 値 と と も に 検 出 す る 場 合 は、 > tkramp(data, start=c(1996,1), frequency=4,tt=c(2008,1,2010,12), ls=c(2001,3))• tcramp

も 同 様 で あ る。

>

tcramp(data,

start=c(1996,1),

quency=4,tt=c(2008,1,2010,12))[2]: zz <- x12siml(shouhi,sorder=2,trend=10,frequency=4,start=c(1994,1),
vat=c(2014,1),vat=c(1997,1),
rp=c(2008,3,2009,1),rp=c(2009,1,2009,4))#zz <- x12siml(shouhi,sorder=2,trend=10,rp=c(2008,2,2009,1),vat=c(2014,1))## VAT ダミーをトレンドに入れたい場合#plot.ts(zz$trend+zz$dummy[,"vat"])[1] "make mat"
vat

rp

rp

11664.552 95997.246 -4679.330

4536.477

vat

vat
vat

rp

rp
17

fre-

1698.4904

731.3393 -1095.7953

1823.0909

Residual Standard Error=753.1881
R-Square=0.1588
F-statistic (df=4, 6)=0.2832
p-value=0.8788
Estimate

Std.Err t-value Pr(>|t|)vat 11664.552 140435.672

0.08310.9365vat 95997.246

94214.217

1.01890.3476rp-4679.329

9227.136 -0.5071

0.6302

rp4536.477

9897.715

0.66280.4583

Residual Standard Error=628.1707
R-Square=0.3265
F-statistic (df=4, 103)=12.4845
p-value=0
Estimate

Std.Err t-value Pr(>|t|)vat1698.4904

274.1562

6.19530.0000vat

731.3393

277.7734

2.63290.0098rp-1095.7953

907.0976 -1.2080

0.2298rp

1823.0909 1101.7570

1.65470.101018

[9]: zz <-␣,→outlier(shouhi,sorder=2,trend=10,type="ao",tt=c(2012,1,2015,1),vat=c(2014,1))## VAT ダミーをトレンド入れたい場合
plot.ts(zz$trend+zz$dummy[,"vat"])[1] "make mat"X
14613.06
X
1708.676
Residual Standard Error=669.3823
R-Square=0.0034
F-statistic (df=1, 9)=0.0307
19

p-value=0.8648
Estimate

Std.Err t-value Pr(>|t|)

X 14613.06 83426.11

0.17520.8648Residual Standard Error=651.6736
R-Square=0.2541
F-statistic (df=1, 106)=36.1055
p-value=0
Estimate

Std.Err t-value Pr(>|t|)

X 1708.676 284.3632

ao

6.00880vat-32341.1 187586.3
aovat

403.5853 1731.1860
aovat-13082.94

83211.63

ao

vat

(中略)
ao

vat

20261.99 78383.67
aovat-64.45557 1711.58885
aovat

22799.13 99500.00
aovat

442.8487 1742.4985

20

21

[3]: source("x12siml92.R")
zz <- x12siml(shouhi,sorder=2,trend=40,frequency=4,start=c(1994,1),
vat=c(2014,1),vat=c(1997,1))
trendvat <- zz$trend + apply(cbind(zz$dummy),1,sum)plot.ts(trendvat)
shouhi2 <- shouhi-apply(cbind(zz$dummy),1,sum)zz2 <- x12siml(shouhi2,sorder=2,trend=40,frequency=4,start=c(1994,1),
rp=c(2008,3,2009,1),rp=c(2009,1,2009,4))[1] "make mat"
vat

vat

6481.018 1414.665
vat

vat
22

1707.0046

724.8809

Residual Standard Error=500.7721
R-Square=0.2517
F-statistic (df=2, 38)=6.3921
p-value=0.004
Estimate

Std.Err t-value Pr(>|t|)vat 6481.018 1844.716

3.51330.0012

vat 1414.665 1714.464

0.82510.4144

Residual Standard Error=639.6484
R-Square=0.2953
F-statistic (df=2, 105)=22.0051
p-value=0
Estimate

Std.Err t-value Pr(>|t|)vat 1707.0046 279.1154

6.11580.0000vat

2.56340.0118724.8809 282.7774

23

rp

rp

-1649.172

1954.881

rprp-1060.731

1836.088

Residual Standard Error=508.1073
R-Square=0.1263
F-statistic (df=2, 38)=2.7471
p-value=0.0769
Estimate
rp -1649.172
rp

Std.Err t-value Pr(>|t|)935.7152 -1.7625

1954.881 1009.1918

1.93710.08600.060224

Residual Standard Error=627.6136
R-Square=0.0373
F-statistic (df=2, 105)=2.0327
p-value=0.1361
Estimate
rp -1060.731
rp

Std.Err t-value Pr(>|t|)906.6988 -1.1699

1836.088 1103.1782

1.66440.24470.099025

[4]: #source("x12siml7m.R")
zz <- x12siml(shouhi,sorder=2,trend=40,frequency=4,start=c(1994,1),
vatt=c(2014,1),#vatt=c(1997,1),# rp=c(2008,3,2009,1),rp=c(2009,1,2009,4))
X
6410.93X
1707.833
Residual Standard Error=498.7189
R-Square=0.2383
F-statistic (df=1, 39)=12.2032
26

p-value=0.0012
EstimateX

Std.Err t-value Pr(>|t|)6410.93 1835.204

3.49330.0012

Residual Standard Error=656.2426
R-Square=0.2513
F-statistic (df=1, 106)=35.5695
p-value=0
Estimate

Std.Err t-value Pr(>|t|)

X 1707.833 286.3562

5.964

0

27

[5]: plot.ts(zz$trend)[11]: juch <- read.csv("juchuuw.csv")[,1]
plot.ts(juch)28[12]: #juch[40] <- juch[40]*1.3
zz <-␣
,→outlier(juch,tt=c(2008,1,2009,10),ilog=T,start=c(2005,4),frequency=12,trend=18,sorder=3,t[1] "make mat"X
-0.1299419
X-0.09596938[1] "-1552.08081898594" "-943.090978293566" "2008"[4] "1""2008""2"X-0.1522739
29

X-0.3235478[1] "-1558.82185151989" "-946.269588459274" "2008"[4] "1""2008""3"X-0.1751104
X-0.2371553[1] "-1566.99742143887" "-943.827425459942" "2008"[4] "1""2008""4"(中略)X
0.06665927
X-0.07261937
X
0.1003971
X-0.04969224
X
0.07624362
X-0.04833362
X
0.03700928
X-0.02554626
X
0.01324179
X-0.002374418
X-0.02469234
X-0.389145
X30

-0.1977416
Residual Standard Error=0.0179
R-Square=0.6918
F-statistic (df=1, 17)=38.1612
p-value=0
Estimate Std.Err t-value Pr(>|t|)X-0.3891

0.063 -6.1775

0Residual Standard Error=0.1186
R-Square=0.0064
F-statistic (df=1, 181)=1.1697
p-value=0.2809
Estimate Std.Err t-value Pr(>|t|)X-0.1977

0.1828 -1.0815

0.2809

31

32

[ ]:
[38]: x <- rnorm(100)x2 <- rnorm(100)y <- x+rnorm(100)lsfit(cbind(x,x2)[1:15,],y[1:15],inter=F)$coef
x[16:100] <- 0
lsfit(cbind(x,x2),y,inter=F)$coefx1.08999511636473 x2-0.216014865209377

x1.09077684498083 x2-0.04949255476582

33

[7]: zz <- x12siml(kakeim[,2],sorder=3,trend=25,frequency=12,start=c(2002,1),
vat=c(2014,3),#vat=c(1997,1),
rp=c(2008,12,2009,4),rp=c(2009,4,2009,12))[1] "make mat"
vat

rp

rp

11988.9364-160.1981

676.9128

vat

rp

rp

1059.58764-30.28929

474.97769

Residual Standard Error=49.5741
R-Square=0.5755
F-statistic (df=3, 22)=9.94
p-value=2e-04
Estimate

Std.Err t-value Pr(>|t|)vat 11988.9364 3103.2146

3.86340.0008rp-160.1981

177.4699 -0.9027

0.3765rp676.9128191.5743

0.00193.5334

Residual Standard Error=360.6754
R-Square=0.1838
F-statistic (df=3, 185)=13.885
p-value=0
Estimate

Std.Err t-value Pr(>|t|)vat 1059.5876

164.5208

6.44040.0000rp-30.2893

770.0244 -0.0393

0.9687rp474.9777 1055.2663

0.45010.653234

[8]: y <- (120:1)/30+rnorm(120)/3
z1 <- x12siml(y,trend=10,sorder=0,start=c(1,1),frequency=1)
y[80] <- y[80]+3#y[31] <- y[31]-3
z2 <- x12siml(y,trend=10,sorder=0,start=c(1,1),frequency=1) #,aot=c(80,1))
plot.ts((z1$trend-z2$trend)[-(80)])[1] "make mat"

35

36

[9]: shouhi3 <- shouhi2shouhi3[99] <- shouhi3[99]*1.1# 季節性をなしにして、Z 全体に回帰させる場合
zz2 <- x12siml(shouhi3,sorder=0,trend=40,frequency=4,start=c(1994,1) ,#rp=c(2008,3,2009,1),rp=c(2009,1,2009,4),
aot=c(2018,3))

shouhi3[80]-shouhi2[80]# 実際の異常値よりも大きい値が推定されてしまう。（１０００ぐらい大きい）37

[1] "make mat"X
9033.126X
6943.179
Residual Standard Error=542.656
R-Square=0.2676
F-statistic (df=1, 39)=14.2473
p-value=5e-04
Estimate

Std.Err t-value Pr(>|t|)

X 9033.126 2393.157

3.7746

5e-04

Residual Standard Error=1318.145
R-Square=0.3592
F-statistic (df=1, 99)=55.4909
p-value=0
Estimate

Std.Err t-value Pr(>|t|)

X 6943.179 932.0679

7.44920038

[10]: shouhi3 <- shouhi2shouhi3[80] <- shouhi3[80]*1.1
zz2 <- x12siml(shouhi3,sorder=2,trend=40,frequency=4,start=c(1994,1) ,#rp=c(2008,3,2009,1),rp=c(2009,1,2009,4),
aot=c(2013,4))# 季節性を仮定して、季節周波数以外で回帰すると、正しい値を推定することができる。
[1] "make mat"X
11454.64
X
7572.638
39

Residual Standard Error=488.4421
R-Square=0.6272
F-statistic (df=1, 39)=65.6134
p-value=0
Estimate

Std.Err t-value Pr(>|t|)

X 11454.64 1414.116

8.10020

Residual Standard Error=636.3842
R-Square=0.7052
F-statistic (df=1, 106)=253.5785
p-value=0
Estimate

Std.Err t-value Pr(>|t|)

X 7572.638 475.5443 15.9241

0

40

[ ]:

41