시계열함수 정리

Summary

본 내용은 코세라 PTSA 과정 중 Lesson 6.3.1 까지 나온 강의 내용을 일부 정리했습니다. 자세한 내용은 해당 강의를 참고하세요.

R 에서는 다양한 시계열 함수 arima, sarima, sarima.for, forecast, auto.arima 가 있다. 이에 대해서 간단한 설명과 사용법을 정리한다. 특히 여기서 사용되는 모든 코드는 코세라 PTSA 강좌의 내용을 일부 수정해서 사용한다.

Modeling 일반적인 과정은 아래와 같다. (출처: https://blog.naver.com/skkong89/223100035178)

• Time plot
• Transformation
• Differencing (seasonal or non-seasonal)
• Ljung-Box test
• ACF => Closers spkies => MA order
• ACF => Spikes around seasonal lags => SMA order
• PACF => Closers spikes => AR order
• PACF => Spikes around seasonal lags => SAR order
• Fit few different models
• Compare AIC, choose a model with minimum AIC
• The parsimony principle
• Time plot, ACF and PACF of residuals
• Ljung-Box test for residuals

데이터 생성 (simulation)

arima.sim(), sarima.sim() 을 이용해서 데이터를 생성할 수 있다.

library(astsa)

set.seed(2023)

n <- 1000 # 데이터 개수

# AR(2) 시뮬레이션
x <- arima.sim(n, model = list(ar = c(0.5, 0.3)), sd=10)
plot(x)

# MA(2) 시뮬레이션
x <- arima.sim(n, model = list(ma = c(0.5, 0.3)), sd=10)
plot(x)

# 예제내용을 그대로 인용함
## AR(2) with mean 50 [n = 500 is default]
y = sarima.sim(ar=c(1.5,-.75)) + 50
tsplot(y) # 시계열에 대한 플롯 생성

## ARIMA(0,1,1) with drift ['mean' refers to the innovations] 
tsplot(sarima.sim(ma=-.8, d=1, mean=.1))

# SARIMA(0, 1, 1, 1, 1, 0)_4 시뮬레이션
tsplot(
  sarima.sim(n = 1000, ar = NULL, d = 1, ma = c(-0.6795), 
             sar = c(-0.322), D = 1, sma = NULL, S = 4)
)

데이터 가시화

• acf: MA(q) 모델에서 지연 q 이후에 0 에 가까워진다.
• pacf: AR(p) 모델은 지연 p 이후에 0 에 가까워진다.

# jj # Johnson & Johnson 분기별 주식 수익
x <- jj

par(mfrow = c(3, 1))
plot(x)
acf(x) # 자기상관함수 Autocorrelation function
pacf(x) # 편자기상관함수 Partial Autocorrelation function

추세를 제거하기 위해 차분을 사용한다. 분산을 제거하기 위해 log 를 사용한다. jj 데이터셋은 시간이 지남에 따라서 추세(평균) 이 증가하고, 분산이 커지는 형태를 볼 수 있다. 또한 계절성을 가지고 있음을 볼 수 있다.

par(mfrow = c(3, 1))
plot(log(jj)) # 분산 제거
plot(diff(log(jj))) # 비계절성 차분, 추세 제거
plot(diff(diff(log(jj)), 4)) # 계절성 차분, 계절성 제거

log 적용하고차분된 데이터에 대해서 acf, pacf 를 확인한다.

data <- diff(diff(log(jj)), 4)

par(mfrow = c(3, 1))
plot(data)
acf(data) # q 는 0 또는 1이 됨을 볼 수 있다. 그리고 대문자 Q 도 0 또는 1이 될거야.
pacf(data) # p 는 0 또는 1, 대문자 P 도 0 또는 1이 될거야.

모델 계수 추정

arima()

용도: 모델의 차수가 주어졌을 때, 모델의 계수를 추정한다.

기본 통계 패키지인 stats 에 포함된다.

# SARIMA(0, 1, 1, 1, 1, 0)_4 

p <- 0; d <- 1; q <- 1 # 비계절성 AR, MA 성분
P <- 1; D <- 1; Q <- 0 # 계절성 AR, MA 성분
span <- 4 # 계절 주기

model <- arima(x = log(jj), order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = span))
model

## 
## Call:
## arima(x = log(jj), order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = span))
## 
## Coefficients:
##           ma1     sar1
##       -0.6796  -0.3220
## s.e.   0.0969   0.1124
## 
## sigma^2 estimated as 0.007913:  log likelihood = 78.46,  aic = -150.91

# 잔차
sse <- sum(model$residuals^2)
sse

## [1] 0.6251634

# 또는 아래와 같은 방식
# sum(resid(model)^2)

## 관측치와 모델치 플롯
plot(log(jj))
lines(fitted(model), col='red')

sarima()

용도: 모델의 차수가 주어졌을 때, 모델의 계수를 추정한다.

• astsa 패키지에 포함된다.
• 내부적으로는 arima() 를 사용한다.
• 추정치에 대한 t.value, p.value 는 model$ttable 을 통해서 확인한다.

위의 내용을 그대로 sarima() 로 코딩한다.

model_s <- sarima(log(jj), p,d,q, P,D,Q, span)

## initial  value -2.237259 
## iter   2 value -2.429075
## iter   3 value -2.446738
## iter   4 value -2.455821
## iter   5 value -2.459761
## iter   6 value -2.462511
## iter   7 value -2.462602
## iter   8 value -2.462749
## iter   9 value -2.462749
## iter   9 value -2.462749
## iter   9 value -2.462749
## final  value -2.462749 
## converged
## initial  value -2.411490 
## iter   2 value -2.412022
## iter   3 value -2.412060
## iter   4 value -2.412062
## iter   4 value -2.412062
## iter   4 value -2.412062
## final  value -2.412062 
## converged

model_s

## $fit
## 
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S), 
##     include.mean = !no.constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
##         REPORT = 1, reltol = tol))
## 
## Coefficients:
##           ma1     sar1
##       -0.6796  -0.3220
## s.e.   0.0969   0.1124
## 
## sigma^2 estimated as 0.007913:  log likelihood = 78.46,  aic = -150.91
## 
## $degrees_of_freedom
## [1] 77
## 
## $ttable
##      Estimate     SE t.value p.value
## ma1   -0.6796 0.0969 -7.0104  0.0000
## sar1  -0.3220 0.1124 -2.8641  0.0054
## 
## $AIC
## [1] -1.910297
## 
## $AICc
## [1] -1.908298
## 
## $BIC
## [1] -1.820318

## XXX: 추정치, t.value, p.value 확인
# sarima 모델에 대해서만 확인이 가능하다.
model_s$ttable

##      Estimate     SE t.value p.value
## ma1   -0.6796 0.0969 -7.0104  0.0000
## sar1  -0.3220 0.1124 -2.8641  0.0054

## 잔차
# 내부적으로 사용되는 arima() 결과가 fit 에 저장되기 때문에 그대로 fit 안의 residuals 를 사용할 수 있다.
sse <- sum(model_s$fit$residuals^2)
sse

## [1] 0.6251634

잔차 분석

Box.test()

잔차내에 AR 성분이 있는지 검정한다. 모델의 잔차에 대해서 AR 성분이 없어야 한다.

## 원본 데이터셋에 대해서 자기상관 여부를 확인
Box.test(jj, lag = log(length(jj)))

## 
##  Box-Pierce test
## 
## data:  jj
## X-squared = 253.48, df = 4.4308, p-value < 2.2e-16

## log(jj)
Box.test(log(jj), lag = log(length(jj)))

## 
##  Box-Pierce test
## 
## data:  log(jj)
## X-squared = 280.52, df = 4.4308, p-value < 2.2e-16

## diff(log(jj))
Box.test(diff(log(jj)), lag = log(length(jj)))

## 
##  Box-Pierce test
## 
## data:  diff(log(jj))
## X-squared = 79.569, df = 4.4308, p-value = 4.441e-16

## diff(diff(log(jj)),4)
Box.test(diff(diff(log(jj)),4), lag = log(length(jj)))

## 
##  Box-Pierce test
## 
## data:  diff(diff(log(jj)), 4)
## X-squared = 20.95, df = 4.4308, p-value = 0.0004938

## 잔차내 AR 성분 검정
Box.test(model$residuals, lag = log(length(jj)))

## 
##  Box-Pierce test
## 
## data:  model$residuals
## X-squared = 2.4914, df = 4.4308, p-value = 0.7079

모델 예상(예측) 수행

forecast() 또는 predict()

용도: arima 모델을 이용해서 예측한다.

• forecast 패키지에 포함된다.

library(forecast)

plot(forecast(model))

forecast(model)

##         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1981 Q1       2.910254 2.796250 3.024258 2.735900 3.084608
## 1981 Q2       2.817218 2.697507 2.936929 2.634135 3.000300
## 1981 Q3       2.920738 2.795579 3.045896 2.729325 3.112151
## 1981 Q4       2.574797 2.444419 2.705175 2.375401 2.774194
## 1982 Q1       3.041247 2.868176 3.214317 2.776558 3.305935
## 1982 Q2       2.946224 2.762623 3.129824 2.665431 3.227016
## 1982 Q3       3.044757 2.851198 3.238316 2.748734 3.340780
## 1982 Q4       2.706534 2.503505 2.909564 2.396028 3.017041

predict(model, 8)

## $pred
##          Qtr1     Qtr2     Qtr3     Qtr4
## 1981 2.910254 2.817218 2.920738 2.574797
## 1982 3.041247 2.946224 3.044757 2.706534
## 
## $se
##            Qtr1       Qtr2       Qtr3       Qtr4
## 1981 0.08895758 0.09341102 0.09766159 0.10173473
## 1982 0.13504742 0.14326433 0.15103487 0.15842473

sarima.for()

sarima 모델에 대해서 예측을 수행한다.

sarima.for(log(jj), n.ahead = 8, p,d,q, P,D,Q, span)

## $pred
##          Qtr1     Qtr2     Qtr3     Qtr4
## 1981 2.910254 2.817218 2.920738 2.574797
## 1982 3.041247 2.946224 3.044757 2.706534
## 
## $se
##            Qtr1       Qtr2       Qtr3       Qtr4
## 1981 0.08895758 0.09341102 0.09766159 0.10173473
## 1982 0.13504742 0.14326433 0.15103487 0.15842473

# XXX: 아래처럼 사용할 수 없다.
# plot(forecast(model_s$fit))

모델의 차수 선택

데이터셋 jj 에 대한 차수를 선택한다. 차수를 선택하는 기준으로 가장 작은 AIC 를 사용한다.

## 모델 선택
# ------------------

# AIC 기준으로 모델의 차수를 계산한다.
my_order_model <- function (d, DD, per, data_ts) 
{
  # d <- 0 # 비계절성 차분
  # DD <- 1 # 계절성 차분
  # per <- 7 # 계절성 주기
  
  list_aic <- NULL
  list_sse <- NULL
  list_pvalue <- NULL 

  idx = 1
  
  # 차수 결정
  for(p in 1:2){
    for(q in 1:2){
      for(i in 1:2){
        for(j in 1:2){
          if(p+d+q+i+DD+j<=10){

            model<-arima(x=data_ts, order = c((p-1),d,(q-1)), seasonal = list(order=c((i-1),DD,(j-1)), period=per))
            pval<-Box.test(model$residuals, lag=log(length(model$residuals)))
            sse<-sum(model$residuals^2)

            list_aic[idx] <- model$aic
            list_sse[idx] <- sse
            list_pvalue[idx] <- pval$p.value

            cat('idx:', idx, ', (', p-1,d,q-1,i-1,DD,j-1, ')_', per, 'AIC=', model$aic, ' SSE=',sse,' p-VALUE=', pval$p.value,'\n')

            idx = idx + 1

          }
        }
      }
    }
  }
  
  return(list(aic = list_aic, sse = list_sse, pvalue = list_pvalue))
}


model_order_result <- my_order_model(1, 1, 4, log(jj))

## idx: 1 , ( 0 1 0 0 1 0 )_ 4 AIC= -124.0685  SSE= 0.9377871  p-VALUE= 0.0002610795 
## idx: 2 , ( 0 1 0 0 1 1 )_ 4 AIC= -126.3493  SSE= 0.8856994  p-VALUE= 0.0001606542 
## idx: 3 , ( 0 1 0 1 1 0 )_ 4 AIC= -125.9198  SSE= 0.8908544  p-VALUE= 0.0001978052 
## idx: 4 , ( 0 1 0 1 1 1 )_ 4 AIC= -124.3648  SSE= 0.8854554  p-VALUE= 0.000157403 
## idx: 5 , ( 0 1 1 0 1 0 )_ 4 AIC= -145.5139  SSE= 0.6891988  p-VALUE= 0.03543717 
## idx: 6 , ( 0 1 1 0 1 1 )_ 4 AIC= -150.7528  SSE= 0.6265214  p-VALUE= 0.6089542 
## idx: 7 , ( 0 1 1 1 1 0 )_ 4 AIC= -150.9134  SSE= 0.6251634  p-VALUE= 0.707918 
## idx: 8 , ( 0 1 1 1 1 1 )_ 4 AIC= -149.1317  SSE= 0.6232876  p-VALUE= 0.6780876 
## idx: 9 , ( 1 1 0 0 1 0 )_ 4 AIC= -139.8248  SSE= 0.7467494  p-VALUE= 0.03503386 
## idx: 10 , ( 1 1 0 0 1 1 )_ 4 AIC= -146.0191  SSE= 0.6692691  p-VALUE= 0.5400205 
## idx: 11 , ( 1 1 0 1 1 0 )_ 4 AIC= -146.0319  SSE= 0.6689661  p-VALUE= 0.5612964 
## idx: 12 , ( 1 1 0 1 1 1 )_ 4 AIC= -144.3766  SSE= 0.6658382  p-VALUE= 0.5459445 
## idx: 13 , ( 1 1 1 0 1 0 )_ 4 AIC= -145.8284  SSE= 0.667109  p-VALUE= 0.2200484 
## idx: 14 , ( 1 1 1 0 1 1 )_ 4 AIC= -148.7706  SSE= 0.6263677  p-VALUE= 0.594822 
## idx: 15 , ( 1 1 1 1 1 0 )_ 4 AIC= -148.9175  SSE= 0.6251104  p-VALUE= 0.7195469 
## idx: 16 , ( 1 1 1 1 1 1 )_ 4 AIC= -144.4483  SSE= 0.6097742  p-VALUE= 0.3002702

list_aic <- model_order_result$aic 
list_sse <- model_order_result$sse
list_pvalue <- model_order_result$pvalue

# 모델 선택을 위해서, AIC, SSE 값을 확인한다.
plot(list_aic, type='l')

# plot(list_sse, type='l')


## 가장 작은 AIC 를 가진 차수 선택
#idx: 6 , ( 0 1 1 0 1 1 )_ 4 AIC= -150.7528  SSE= 0.6265214  p-VALUE= 0.6089542 
#idx: 7 , ( 0 1 1 1 1 0 )_ 4 AIC= -150.9134  SSE= 0.6251634  p-VALUE= 0.7079173  제일 작다.
#idx: 8 , ( 0 1 1 1 1 1 )_ 4 AIC= -149.1317  SSE= 0.6232876  p-VALUE= 0.6780876 

model_s <- sarima(log(jj), 0,1,1, 1,1,0, 4)

## initial  value -2.237259 
## iter   2 value -2.429075
## iter   3 value -2.446738
## iter   4 value -2.455821
## iter   5 value -2.459761
## iter   6 value -2.462511
## iter   7 value -2.462602
## iter   8 value -2.462749
## iter   9 value -2.462749
## iter   9 value -2.462749
## iter   9 value -2.462749
## final  value -2.462749 
## converged
## initial  value -2.411490 
## iter   2 value -2.412022
## iter   3 value -2.412060
## iter   4 value -2.412062
## iter   4 value -2.412062
## iter   4 value -2.412062
## final  value -2.412062 
## converged

y_pred <- sarima.for(log(jj), n.ahead = 8, 0,1,1, 1,1,0, 4)

## 원본 데이터셋에 예측치를 추가하고, 시계열로 만든 후 플롯한다.
t1 <- c(jj, exp(y_pred$pred))
t1 <- ts(t1, start = c(1960), frequency = 4)

plot(t1) # 원본 데이터셋 + 예측치
lines(jj, col='blue') # 원본 데이터셋

auto.arima() 사용

위에서 차수를 선택하고, 해당 차수에 대한 모델 계수를 구하는 과정을 거쳤다. 이것을 하나의 함수로 구현한 것이 auto.arima() 이다.

model <- auto.arima(log(jj), ic='aic', approximation = FALSE)
model

## Series: log(jj) 
## ARIMA(2,0,2)(0,1,0)[4] with drift 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2   drift
##       0.6797  -0.6133  -0.4031  0.7998  0.0386
## s.e.  0.1650   0.1693   0.1228  0.1270  0.0035
## 
## sigma^2 = 0.007573:  log likelihood = 83.92
## AIC=-155.83   AICc=-154.68   BIC=-141.54

plot(forecast(model, h=8)) # 관측치 + 예측치 표시

# fitted(model)
# resid(model)
# forecast(model)

plot(log(jj)) # 관측치 표시
lines(fitted(model), col='red') # 예측치 표시

# predict(model, n.ahead = 8) # 사용 불가

HoltWinters

## 모델 피팅
AirPassengers.hw <- HoltWinters(log10(AirPassengers))
AirPassengers.hw

## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = log10(AirPassengers))
## 
## Smoothing parameters:
##  alpha: 0.326612
##  beta : 0.005744246
##  gamma: 0.8207255
## 
## Coefficients:
##             [,1]
## a    2.680598830
## b    0.003900787
## s1  -0.031790733
## s2  -0.061224237
## s3  -0.015941495
## s4   0.006307818
## s5   0.014138008
## s6   0.067260071
## s7   0.127820295
## s8   0.119893006
## s9   0.038321663
## s10 -0.014181699
## s11 -0.085995400
## s12 -0.044672707

## 계수 확인
AirPassengers.hw$coefficients

##            a            b           s1           s2           s3           s4 
##  2.680598830  0.003900787 -0.031790733 -0.061224237 -0.015941495  0.006307818 
##           s5           s6           s7           s8           s9          s10 
##  0.014138008  0.067260071  0.127820295  0.119893006  0.038321663 -0.014181699 
##          s11          s12 
## -0.085995400 -0.044672707

## 예측
# predict(AirPassengers.hw, n.ahead = 1)
forecast(AirPassengers.hw, h = 1)

##          Point Forecast    Lo 80   Hi 80    Lo 95    Hi 95
## Jan 1961       2.652709 2.630898 2.67452 2.619351 2.686066

## 예측
AirPassengers.forecast <- forecast(AirPassengers.hw)
plot(AirPassengers.forecast, xlim=c(1949, 1963))

# plot(exp(AirPassengers.forecast$x))
# lines(exp(AirPassengers.forecast$fitted), col='red')

기타

• ts: 시계열 객체 생성
• frequency: 시계열 객체의 주기
• start, end: 시계열 객체의 시작과 끝 시간
• window: 시간 범위
• stl: LOESS에 의한 시계열 객체 분해
• stats::lag: 시계열 객체의 시작 시간을 이동함

## AirPassengers 가 'ts' 객체임을 확인
class(AirPassengers)

## [1] "ts"

class(EuStockMarkets)

## [1] "mts"    "ts"     "matrix" "array"

## 주기 확인
frequency(AirPassengers)

## [1] 12

## 시작, 끝 시간 
start(AirPassengers)

## [1] 1949    1

end(AirPassengers)

## [1] 1960   12

## 시간 범위: 1955년 2월 ~ 1960년 11월
x <- window(AirPassengers, start = c(1955, 2), end = c(1960, 11))
x

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1955     233 267 269 270 315 364 347 312 274 237 278
## 1956 284 277 317 313 318 374 413 405 355 306 271 306
## 1957 315 301 356 348 355 422 465 467 404 347 305 336
## 1958 340 318 362 348 363 435 491 505 404 359 310 337
## 1959 360 342 406 396 420 472 548 559 463 407 362 405
## 1960 417 391 419 461 472 535 622 606 508 461 390

plot(x)

## stl
plot(stl(AirPassengers, 'periodic'))

## 관심 대상은, 측정치 자체가 아니라 다음 측정치로의 변화정도이다. 히스토그램을 통해 확인한다.
hist(AirPassengers)

hist(diff(AirPassengers))

a <- 1:10
b <- 1:10
a.ts <- ts(a)
b.ts <- ts(b)
stats::lag(a.ts, k = 1) # 이전 시간으로 시작 시간을 옮긴다.

## Time Series:
## Start = 0 
## End = 9 
## Frequency = 1 
##  [1]  1  2  3  4  5  6  7  8  9 10

library(data.table) # for shift()

x <- sample(1:10, 10) # 샘플 데이터 생성
x

##  [1]  2  1  7 10  9  3  8  6  4  5

dplyr::lag(x, 1) # 시차 1로 지연

##  [1] NA  2  1  7 10  9  3  8  6  4

shift(x, 1) # 시차 1로 지연

##  [1] NA  2  1  7 10  9  3  8  6  4

diff(x, 1) # 현재 시점에서 바로 이전 시점을 뺀다.

## [1] -1  6  3 -1 -6  5 -2 -2  1

x - shift(x, 1) # 현재 시점에서 바로 이전 시점을 뺀다.

##  [1] NA -1  6  3 -1 -6  5 -2 -2  1

stats::lag(ts(x), 1) # 데이터 포인트의 시작 날짜/시간 정보를 이동한다.

## Time Series:
## Start = 0 
## End = 9 
## Frequency = 1 
##  [1]  2  1  7 10  9  3  8  6  4  5

stats::lag(ts(x), 1) - ts(x)

## Time Series:
## Start = 1 
## End = 9 
## Frequency = 1 
## [1] -1  6  3 -1 -6  5 -2 -2  1

롤링 윈도

주기성을 제거한다.

[see 닐슨 2021:124]

x <- rnorm(n = 100) + 1:100
plot(x, type='l')
lines(filter(x, rep(1/3, 3)), col='red', lty=2)
lines(filter(x, rep(1/10, 10)), col='blue', lty=3)

자기상관함수 ACF

a <- AirPassengers

print(acf(a))

## 
## Autocorrelations of series 'a', by lag
## 
## 0.0000 0.0833 0.1667 0.2500 0.3333 0.4167 0.5000 0.5833 0.6667 0.7500 0.8333 
##  1.000  0.948  0.876  0.807  0.753  0.714  0.682  0.663  0.656  0.671  0.703 
## 0.9167 1.0000 1.0833 1.1667 1.2500 1.3333 1.4167 1.5000 1.5833 1.6667 1.7500 
##  0.743  0.760  0.713  0.646  0.586  0.538  0.500  0.469  0.450  0.442  0.457

cor(a, shift(a, 1), use = 'pairwise.complete.obs')

## [1] 0.9601946

x <- rnorm(100)
x <- ts(x)
plot(x)

acf(x)

pacf(x)

auto.arima(x)

## Series: x 
## ARIMA(0,0,0) with zero mean 
## 
## sigma^2 = 0.9454:  log likelihood = -139.09
## AIC=280.17   AICc=280.22   BIC=282.78

x <- 1:100
plot(x)

acf(x, lag.max = 30)

pacf(x)

2차원 시각화

[see 닐슨 2021:146]

AirPassengers 데이터를 사용해서 계절성과 추세를 확인

a <- t(matrix(AirPassengers, nrow=12, ncol=12)) # 좌우 변환


## 2차원 시각화: 아래 2개는 동일
matplot(matrix(AirPassengers, nrow=12, ncol=12), type='l', lty=1, lwd=2.5, xaxt = 'n')

library(forecast)
seasonplot(AirPassengers)

## 연도별 월별 곡선
matplot(t(matrix(AirPassengers, nrow=12, ncol=12)), type='l', lty=1, lwd=2.5, xaxt = 'n')

monthplot(AirPassengers)

3차원 시각화

#install.packages('plotly')
library(plotly)

## 필요한 패키지를 로딩중입니다: ggplot2

## 
## 다음의 패키지를 부착합니다: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library(data.table)

months <- 1:12
ap <- data.table(matrix(AirPassengers, nrow = 12, ncol = 12))
names(ap) <- as.character(1949:1960)
ap[, month := months]
ap <- melt(ap, id.vars = 'month')
names(ap) <- c('month', 'year', 'count')

p <- plot_ly(ap, x = ~month, y = ~year, z = ~count, color = ~as.factor(month)) %>%
  add_markers() %>%
  layout(scene = list(xaxis = list(title = 'Month'),
                      yaxis = list(title = 'Year'),
                      zaxis = list(title = 'PassengerCount')))

p

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## 랜덤워크
x <- rnorm(100); x <- cumsum(x)
y <- rnorm(100); y <- cumsum(y)
z <- rnorm(100); z <- cumsum(z)

temp <- data.frame(x = x, y = y, z = z)
p <- plot_ly(temp, x= ~x, y= ~y, z= ~z, type='scatter3d', mode='lines',
             line=list(color='#1f77b4', width=1))

p

References

https://www.coursera.org/learn/practical-time-series-analysis/home/welcome

닐슨. 2021. Practical time series analysis. 한빛미디어.