x12simldoc92
2023 年 2 月 2 日

1

Seasonal Adjustment by Using X12SIML

Author : Seisho Sato (University of Tokyo)
Assisted by Naoto Kunitomo (Institute of Statistical Mathematics)
February 2, 2023

1.1

X12SIML (Ver.92)• Seasonal adjustment based on the SIML method.
• References : Sato and Kunitomo (2020, 2021) :http://www.mims.meiji.ac.jp/publications/2020-ds/SDS-15.pdf
http://www.mims.meiji.ac.jp/publications/2021-ds/SDS-20.pdf• It is like X-12ARIMA of the U.S. Census Bureau with outlier detections.
• Works in R.
• I has auto-outlier detections such as AO, LS, RAMP, and others.
• It is possible to define regression variables in the program.
• Install source program (“x12siml8.R”). No need for other libraries.

[1]: source("x12siml92.R")

1.2

How to use x12siml“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, …)• In reg, you can insert user defined explanatory variables.
• In trend and sorder, insert positive integers. The larger values mean more flexibility on
trend and seasonality. As a suggestion, use the number of 10% of data for trend and 1～
3 for sorder.
1

• pb is the cylce number of naive predictions and pa is the cylce number of naive backcasts.
pb and pa can be used to stabilize the estimates of seasonal states when the number ofdata is not large.
• mtype is the SIML filter-type. It is 1 for forward-filter and 2 for the backward-filter.
• In … , you can use the outlier variables in the following list.
• Structural Breaks and Outliers– AO(additive outlier)– LS (level shift)– TC (temporaly change)– RAMP (slope change shift, ramp)– VAT– Holiday– Leap Year
– TCRAMP(RAMP+TC)– TCRAMP1(RAMP+(one blank)+TC)• Outlier variables• AO > x12siml(data,start=c(1993,1),frequency=4,ao=c(2001,2))(AO is setted at 2001Q2.)• LS > x12siml(data,start=c(1993,1),frequency=4,ls=c(2001,2))(LS is setted at 2001Q2.)• RAMP > x12siml(data,start=c(1993,1),frequency=4,rp=c(2001,2,2001,4))(RAMP is setted from 2001Q2 to Q4.)
• TLS > x12siml(data,start=c(1993,1),frequency=4,tls=c(2001,2,2001,4))(TLS(temporary level shift) is setted from 2001Q2 to Q4.)
• TC >

x12siml(data,start=c(1993,1),frequency=4,tc=c(2001,2))(TC is setted at 2001Q2.)• VAT >

x12siml(data,start=c(1993,1),frequency=4,vat=c(2001,1))(VAT dummy (1 at the selected pwriod and -1 in the next period).)• VATT >x12siml(data,start=c(1993,1),frequency=4,vatt=c(2001,1))(VATT dummy (1 at the selected pwriod and -1 in the next period with a trend).)• Holiday >x12siml(data,start=c(1993,1),frequency=4,hol=c(2001,3,2003,4))(Dummies as AO with several periods (the resulting estimates are in the estimated seasonal).)• LeapYear >x12siml(data,start=c(1993,1),frequency=4, ly=c(1996,1))(Set AO at each 4*frequency after the selected period.)• SLS >

x12siml(data,start=c(1993,1),frequency=4, sls=c(2003,4))(Level shift in seasonals at the selected period.)• TCRAMP > x12siml(data,start=c(1993,1),frequency=4,tcrp=c(2001,2,2001,4,2003,3))(Set TCRAMP from 2001Q2 to 2002Q3 (RAMP from 2001Q2 to 2001Q4 and TC from2

2001Q4 to 2003Q3).)• TCRAMP1 > x12siml(data,start=c(1993,1),frequency=4,tcrp1=c(2001,2,2001,4,2003,3))(Set TCRAMP1 from 2001Q2 to 2002Q3 (RAMP from 2001Q2 to 2001Q4 and TC from
2002Q1 to 2003Q3).)[ ]:[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)

3

tt <- sum((reg2[[i]] - start) * c(frequency,␣
,→

1)
) + 1
z[tt:n] <- tcrate^(0:(n - tt))
reg <- cbind(reg, z)},"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)4

},"ly"= {z <- rep(0, n)if(length(reg2[[i]])==1)
tt <- reg2[[i]]+(0:n)*(frequency*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)},

5

"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
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)] <- -1

6z[(tt1 + 1):(tt2 - 1)] <- (((tt1 + 1):(tt2 -␣,→

1)) - tt1)/(tt2 - tt1)
z[tt2] <- 1z[(tt2 + 1):tt3] <- tcrate^(0:(tt3 - tt2 -␣,→1))
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

Seasonal Frequency• We set the seasonal frequency as s (that is, s = 12 for monthly data, s = 4 for qurterlydata）.
• We usually observe seasonal fluctuations at fs :
1 2
[s/2]
fs = ( , , · · · ,)s ss
f12 = (1 2
6, ,··· , )
12 12
12

1 2 3
f7 = ( , , )
7 7 7
1 2
f4 = ( , )4 411• We set the seasonal band with sorder (flexible seasonals) :inds = (2 ∗ fs ∗ n − sorder) : (2 ∗ fs ∗ n + sorder)(≤ n)
ind4 = (n/2 − sorder) : (n/2 + sorder), (n − sorder) : n
We assume that the number 1～n in the x-axis corresponds to 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
12

plot.ts(sfreq,type="h")abline(h=0)

1.4Regression (regsiml)• A) We set the Kn transformed data as Z, the Kn transformed regression variables as
Z (d) .

z = P C −1 y, Zd = P C −1 W• B) From Z, select the seasonal band and construct Zs .
z s = Js z
Js is the selection matrix to find the seasonal bands.
13• C) From Z, select bands except the seasonal bands and construct Z−s .
z−s = J−s z
J−s is the selection matrix chooshing bands except the seasonal bands.
• D) From Z−s , form Zt by selecting the first m terms (m is trend).
zt = Jm z
Jm is the selection transformation of the first m terms.
• E) For Z (d) , construct Zs ,Z−s ,Zt
(d)

(d)

(d)

in the same way.

Zds = Js Zd, Zd−s = J−s Zd, Zdt = Jm Zd• F) Regression A : Regress Zt on Zt(d)with 0 intercept.

Ct = (Zd0t Zdt )−1 Zd0t zt• G) Regression B : Regress Z−s on Z−s with 0 intercept.
(d)C−s = (Zd0−s Zd−s )−1 Zd0−s z−s• H) When the dummy variables include one of AO, VAT, HOL, LY, then, compare theF values to select the results of either Regression A or Regression B. Otherwise,
Regression A is chosen.
•I) By running Regression, obtain the regressison coefficients ”coef” and the residuals”res”.

0
0
C = Ct or C−s , res = Jm
zt − Zdt Ct or J−s
z−s − Zd−s C−s• J) Zs∗ is the residuals of Zs(d)minus coef times Zs .
zs ∗ = zs − Zds C• K) Zt∗ is composed of the first m terms of res and set 0 otherwise,zt ∗ = Jm res• L) The inverse transformations of Zs∗ and Zt∗ are the seasonal Ｓ and trend Tr .
0
S = CP 0 Js0 zs ∗, T r = CP 0 Jm
zt ∗• M) AIC of Regression A is AIC.a and AIC of Regression B is AIC.b . (We multiple n1
instead of n to Log of estimated variances.)14[8]: ## regsiml
mat <- siml.mat(n,m1,type=mtype)## type=1 (前向き), 2 (後ろ向き)
n1 <- dim(mat$K)[1]## 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:15

1. siml.mat(n, m1, type = mtype)2. makeMat21(n, m, type = type)

1.5

Result (x12siml)- In x12siml, we use the results of regsiml to make figures.
Among dummy variables, except “ao”,“vat”,“hol”,“ly” the effects are included in the estimatedtrend. The variables “hol” and “ly”

are included in the estimated seasonals The variables“ao” and “vat’ are included in the estimated noise.’[ ]: 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]}}16

z.adj <- data - zz.seasonalzz.noise <- data-(zz.seasonal+z.trend)

1.6

Outlier Auto-detection
outlier(data, start=c(1998,1), frequency=12,type=“ao”,tt=c(1999,1,2009,12))• We search AO from 1999, January to 2009 December. Search one period to improve AICwith certain value (the default is 4).
• To

detectmultiple

periods,

set

:

>

outlier(data,

start=c(1998,1),

fre-

quency=12,type=“ao”,tt=c(2001,4,2009,12), ao=c(2001,3))– (Some care should be made not to take double counts of period.)• Fo type, you can use “ao”, “ls”, “tc”, “rp” .
• Auto-search of TKRAMP and/or TCRAMP• Auto-search tkramp during 2008, January～2010, December. (It may take time if thesearch length is long.)• TKRAMP is the double RAMP（Takaoka-RAMP, RAMP+RAMP）• Search one period to improve AIC with certain value (the default is 4).
• To use detection with other outliers, use :> tkramp(data, start=c(1996,1), fre-quency=4,tt=c(2008,1,2010,12), ls=c(2001,3))
• The

same

for

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

If you want to include VAT dummy variable to trend.

%(VAT ダミーをトレンドに入れたい場合)#plot.ts(zz$trend+zz$dummy[,"vat"])[1] "make mat"
vat

vat

rp

rp

11664.552 95997.246 -4679.330

4536.477
17

fre-

vat
1698.4904

vat

rp

rp

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# Without seasonals, regress all Z%(季節性をなしにして、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# The estimated value may be larger than the actual outlier.
(abouut 1000.)%(実際の異常値よりも大きい値が推定されてしまう。（１０００ぐらい大きい）)[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))#

Assuming seasonals and regress except the seasonalfrequency bonds, we can estimate the true value.
%(季節性を仮定して、季節周波数以外で回帰すると、正しい値を推定することができる。)[1] "make mat"X
39

11454.64
X
7572.638
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