ARMA 1 (AIC) ARMA 2 ARMA dy(t) = f(t, y), y(0) = y 0 dt f(t, y) y 0 f(t, y) LCR (Fig. 1) k c m t y(t) F = ma F = θ(t) ky(t) c dy(t), ma = m d2 y(t) dt dt 2 m d2 y(t) dt 2 + c dy(t) dy(0) + ky(t) = v(t), y(0) = y 0, = y 1 dt dt v(t) y 0 Y 1 m v(t) v(t) y(t) y(t) y[n] y[n] 2y[n 1] + y[n 2] y[n] y[n 2] m 2 + c + ky[n 1] = v[n] 2 Delta y[n] = a 1 y[n 1] + a 2 y[n 2] + b 0 v[n] (n = 0, 1, 2,, N 1) ( ) /( 2m m a 1 = 2 k 2 + c ) 2 ( a 2 = m 2 + c ) /( m 2 2 + c ) 2 /( m b 0 = 1 2 + c ) 2
y[n] y[n 1], y[n 2] v[n] R k c m y(t) Figure 1: series.r series <- function(n,delta,k,c,m,y0,y1) { a1 <- (2* m/delta^2 - k) / (m/delta^2 + 0.5*c/delta); a2 <- (- m/delta^2 + 0.5*c/delta) / (m/delta^2 + 0.5*c/delta); b0 <- 1 / (m/delta^2 + 0.5*c/delta); y <- numeric(n+2); v <- numeric(n); # for simulation 1 #v <- rnorm(n); # for simulation 2 y[1] <- y0; y[2] <- y[1] + y1*delta; list = 1:n; for(i in list) { y[i+2] <- a1*y[i+1] + a2*y[i] + v[i] return(y); Fig. 2 N = 256, = 1, k = 0.2, c = 0.06, m = 1, y0 = 1, y1 = 0 v[n] = 0 (n = 0, 1,, N 1)
(1) > source("series.r") > y <- series(256,1,0.2,0.06,1,1,0) > plot(y, type="l") > dev.copy2eps(file="simulation1.eps", width=10, height=6) y 0.5 0.0 0.5 1.0 Figure 2: (1) (v[n] = 0) y 0 = y 1 = 0 Fig. 3 (1) > source("series.r") > y <- series(256,1,0.2,0.06,1,0,0) > plot(y, type="l") > dev.copy2eps(file="simulation2.eps", width=10, height=6 )
y 10 0 10 20 Figure 3: (2) 3 ARMA y[n] v[n] y[n] = a i y[n i] + v[n] b i v[n i] ARMA (AutoRegressive Moving Average model) m, {a i m (AR coefficients) l, {b i l (MA coefficients) v[n] 0 σ 2 y[n i] v[n] E{v[n] = 0, E{v[n]v[m] = σ 2 δ mn, E{v[n]y[m] = 0 (n > m) δ mn { 1 (n = m) δ mn = 0 Otherwise ARMA l = 0 y[n] = a i y[n i] + v[n] m (AutoRegressive Model) AR m = 0 (Moving Average Model) MA y[n] = v[n] b i v[n i]
(a) 2 AR y[n] = 0.9 3y[n 1] 0.81y[n 2] + v[n] (b) 2 MA y[n] = v[n] 0.9 2v[n 1] + 0.81v[n 2] (c) (2, 2) ARMA y[n] = 0.9 3y[n 1] 0.81y[n 2] + v[n] 0.9 2v[n 1] + 0.81v[n 2] R > n <- 256 > m <- 2 > a <- c(0.9*sqrt(3), -0.81) > l <- 2 > b <- c(0.9*sqrt(2), -0.81) >sigma2 <- 1 # AR(2) model > source("armodel.r") > y <- armodel(n,m,a,sigma2) > plot(y, type="l") > dev.copy2eps(file="armodel.eps", width=10, height=6) # MA(2) model > source("mamodel.r") > y <- mamodel(n,l,b,sigma2) > plot(y, type="l") > dev.copy2eps(file="mamodel.eps", width=10, height=6) # ARMA(2,2) model > source("armamodel.r") > y <- armamodel(n,m,l,a,b,sigma2) > plot(y, type="l") > dev.copy2eps(file="armamodel.eps", width=10, height=6)
ARmodel.R # AR model armodel <- function(n,m,a,sigma2) { y <- numeric(n+m); v <- rnorm(n, mean=0, sd=sqrt(sigma2)); listi = 1:n; listj = 1:m; for(i in listi) { y[i+m] <- v[i]; for(j in listj) { y[i+m] <- y[i+m] + a[j] * y[i+m-j]; y <- c(y)[m+1:n]; return(y); y 5 0 5 10 Figure 4: (AR model)
MAmodel.R # MA model mamodel <- function(n,l,b,sigma2) { y <- numeric(n); v <- c(numeric(l),rnorm(n, mean=0, sd=sqrt(sigma2))); listi = 1:n; listj = 1:l; for(i in listi) { y[i] <- v[i+l]; for(j in listj) { y[i] <- y[i] - b[j] * v[i+l-j]; return(y); y 4 2 0 2 4 6 Figure 5: (MA model)
ARMAmodel.R # ARMA model armamodel <- function(n,m,l,a,b,sigma2) { y <- numeric(n+m); v <- c(numeric(l),rnorm(n, mean=0, sd = sqrt(sigma2))); listi = 1:n; listj = 1:l; listk = 1:m; for(i in listi) { y[i+m] <- v[i+l]; for(j in listj) { y[i+m] <- y[i+m] - b[j] * v[i+l-j]; for (k in listk) { y[i+m] <- y[i+m] + a[k] * y[i+m-k]; y <- c(y)[m+1:n]; return(y); y 4 2 0 2 4 Figure 6: (ARMA model)
4 ARMA ARMA z 1 D ARMA ( ) ( ) 1 a i z i y[n] = 1 b i z i v[n] AR MA ( ) a(z 1 ) := 1 a i z i ARMA a(z 1 )y[n] = b(z 1 )v[n] ( ) b(z 1 ) := 1 b i z i z 1 g(z 1 ) = { a(z 1 ) 1 b(z 1 ) = g i z i i=0 y[n] = g(z 1 )v[n] = g i v[n i] i=0 ARMA MA g i ARMA i g 0 = 1, g i = a j g i j b i (i > 1) j=1 AR MA a j = 0 (j > m) b i = 0 (i > l) ARMA 5 ARMA ARMA y[n k] E(y[n]y[n k]) = a i E(y[n i]y[n k]) + E(v[n]y[n k]) b i E(v[n i]y[n k]) ARMA MA { 0 (n > m) E(v[n]y[m]) = g i E(v[n]v[m i]) = σ 2 g m n (n m) i=0
ARMA C 0 = C k = ) a i C i + σ (1 2 b i g i a i C k i σ 2 b i g i k (k = 1, 2, ) ARMA AR MA {m, a 1, a 2,, a m, {l, b 1, b 2,, b l σ 2 {g i,2, {C k k=0,1,2, AR C 0 = C k = a i C i + σ 2 a i C k i (k = 1, 2, ) Yule- Walker AR 6 ARMA ARMA C k (DTFT) p(f) = DTFT{C k = k= C k exp( 2πjkf) =< C k, exp(2πjkf) > ( 1 2 f 1 ) 2 C k (power spectral density function) (Spectrum) p(f) = = = = k= C k exp( 2πjkf) E(y[n]y[n k]) exp( 2πjkf) k= {( ) ( ) E g r v[n r] g s v[n k s] exp( 2πjkf) k= r=0 s=0 g r g s E(v[n r]v[n k s]) exp( 2πjkf) k= r=0 s=0 g r = 0 (r < 0) p(f) = σ 2 g r g r k exp( 2πjkf) k= r=0
= σ 2 r g r exp( 2πjrf)g r k exp{ 2πj(k r)f r=0 k= = σ 2 g r exp( 2πjrf)g s exp( 2πjsf) = r=0 s=0 σ 2 2 g r exp( 2πjrf) r=0 z = exp(2πjf) g(z 1 ) { 1 { exp( 2πjrf) = 1 a i exp( 2πjif) 1 b i exp( 2πjif) r=0 ARMA 1 l p(f) = σ 2 b i exp( 2πjif) 2 1 m a i exp( 2πjif) 2 ARMA AR {a i m MA {b i l σ 2 AR(1),AR(2),MA(1) ARMA(1,1) AR MA ARMA R plot(y, type= l ) y acf(y) y spec.pgram(y, type= h ) y AR(2) MA(2) ARMA(2,2)
y 5 5 Series y ACF 0.5 0.5 0 5 10 15 20 Lag Series: y Raw Periodogram spectrum 1e 03 0.0 0.1 0.2 0.3 0.4 0.5 frequency bandwidth = 0.00113 Figure 7: AR model y 4 0 2 4 Series y ACF 0.5 0.0 0.5 1.0 0 5 10 15 20 Lag Series: y Raw Periodogram spectrum 1e 03 1e 01 1e+01 0.0 0.1 0.2 0.3 0.4 0.5 frequency bandwidth = 0.00113 Figure 8: MA model
y 4 2 0 2 Series y ACF 0.2 0.2 0.6 1.0 0 5 10 15 20 Lag spectrum 1e 03 1e 01 1e+01 Series: y Raw Periodogram 0.0 0.1 0.2 0.3 0.4 0.5 frequency bandwidth = 0.00113 Figure 9: ARMA model