SetDirectory@NotebookDirectory[]

시계열과 Mathematica

Wolfram Language Conference 2023
Nov 10, 2023
고려대학교
고려대학교 수학과 김영욱

시계열의 기본 개념

시계열이란?

시계열(time series)이란 시간 순서에 따라 연동되어 나타나는 데이터를 말한다.

특히 이산적으로 주어지는 (수)열 형태인 경우가 많다.

예를 들면 간조᛫만조 시의 수위, 태양의 흑점의 개수, 매일 매일의 KOSPI 지수 종가(closing value).

시계열을 가지고 하고 싶은 것

주어진 시계열 데이터에서 유의미한 정보(통계량)을 추출하고 싶다.

특히 미래의 데이터를 예측하려 한다: Forcasting.

구체적으로 두 가지 문제가 있다.

주어진 데이터를 분류한다. (Classification, 검정)

주어진 데이터에서 특정 양을 유추한다. (Regression, 추정)

시계열의 분석 방법들

시계열 데이터를 시간의 함수로 둔채로 분석하는 time-domain 방법과 주파수별로 분해해서 분석하는 frequency-domain 방법이 있다.

Null

위의 그림은 Wikipedia > Frequency domain 에서 가져왔다.

Parametric 방법과 Non-parametric 방법

Parametric 방법: Model을 가지고 있고 그 모델이 유한개의 parameter로 정의되는 경우. 이 확률과정이 구조를 가지고 있다고 생각한다.

SARIMA 등의 형태를 가진다. 모델의 parameter를 통계적으로 유추할 수 있다.

Non-parametric 방법: 모델을 통하지 않고 직접 통계적 특성(공분산, 스펙트럼 등)을 계산하게 된다.

시계열 분석 도구

고전적 통계 분석 도구: Auto-correlation 등.

수학적 도구: Fourier 분석, PCA 등.

Denoising: Filters 등등.

기계학습: Hidden Markov / State Space model, Support vector machine, Neural networks

시간-주파수 분석법: FFT, STFT, Wavelet 변환 등등.

기타 다양한 방법이 있다.

우리가 알아보려는 것

Financial time series와 관련된 기초적 사항

특히 주가 분석에 쓰이는 technical index들

Wavelet 변환

Neural Network

금융 시계열 Financial Time Series

ts=TimeSeries[];
DateListPlot[ts]

ts["May 24, 2009"]

sbux=FinancialData["SBUX","Close",{{2014,5,1},{2014,9,30}}]

DateListPlot@sbux

RegularlySampledQ[sbux]

dates=DayRange[{2014,5,1},{2014,9,30},"BusinessDay"];
values=sbux["Path"][[All,2]];

prices=TimeSeries[values,{{2014,5,1},{2014,9,30},"BusinessDay"}];
RegularlySampledQ[prices]

prices2=TimeSeries[sbux,TemporalRegularity->True];
RegularlySampledQ[prices2]

{sbux["MinimumTimeIncrement"],prices["MinimumTimeIncrement"],prices2["MinimumTimeIncrement"]}

TimeSeries 다루기

TimeSeriesShift
TimeSeriesWindow
TimeSeriesInsert
TimeSeriesRescale
TimeSeriesResample
TimeSeriesAggregate
TimeSeriesMap
TimeSeriesThread
TimeSeriesMapThread
MovingMap and MovingAverage

TimeSeries[{1.2,0.3,4.7},{{2,4,6}}];
TimeSeriesWindow[%,{3,6}];
%//Normal

Missing 데이터

missd={2,1,3,Missing[],2,1,2,Missing[],6,2,5};
missd//Length

ts=TimeSeries[missd,Automatic];

ListLinePlot[ts]

ts["Values"]

ts[[2,1,1]]==ts["Values"]

DeleteMissing@ts["Values"]

m=Mean[DeleteMissing@ts["Values"]]

ts1=TimeSeries[ts["Values"],Automatic,MissingDataMethod->{"Constant",m}]

ts2=TimeSeries[ts["Values"],Automatic,MissingDataMethod->{"Interpolation",InterpolationOrder->1}]

DateListPlot[#]&/@{ts,ts1,ts2}

Mathematica 시계열의 기본

Mathematica는 시계열에서 세 가지 기본 object를 사용한다: TimeSeries, EventSeries, TemporalData.

TimeSeries, EventSeries: 10.0(2014), 10.1(2015)

TemporalData: 9.0(2012), ... , 12.0(2019)

TimeSeries = 연속 시간, EventSeries = 이산 시간, TemporalData = multi path.

TimeSeries는 시간영역이 연속이든 이산적이든 사용할 수 있으며 이산적인 경우에도 보간법으로 간극을 채워 연속으로 만들 수 있다.
values는 array가 될 수 있다.

v={2,1,6,5,7,4};
t={1,2,5,10,12,15};
ts=TimeSeries[v,{t}];
ListLinePlot[ts]

예: TimeSeries의 연산

ts1=TimeSeries[{2,4,6,8},Automatic];
ts2=TimeSeries[{1,3,5,7},Automatic];
sum12=ts1+ts2;
ListLinePlot[{Labeled[ts1,"ts1"],Labeled[ts2,"ts2"],Labeled[sum12,"sum"]},PlotTheme->"OpenMarkers"]

고전적 계산 방법

데이터의 특성을 보여주는 지표를 만든다.

주가(주식의 종가)의 경우에는 여러 가지 지표가 만들어져 있다

FinancialIndicator[All]

TradingChart[{"GOOGL",{"Jul 1, 2022","Jun 31, 2023"}},{"UltimateOscillator","RelativeStrengthIndex",FinancialIndicator["SimpleMovingAverage",50]},PlotLabel->"Google Stock Price — 2022-23"]

Indicator 구성 예

Yahoo finance에서 본 것

First Try

Better One

날자 간격, 휴일 등을 무시하고 단순한 time series로 보고 그린 것

Final 비교

코드

drawIchimoku[fileName_,start_,end_,clPerioD_,blPerioD_,leadSpanBPerioD_,lagSpanPerioD_]:=
Module[{
filenamE=fileName,starT=start,enD=end,clPeriod=clPerioD,blPeriod=blPerioD,leadSpanBPeriod=leadSpanBPerioD,lagSpanPeriod=lagSpanPerioD,
data,raw,dataDates,dataClose,dataHigh,dataLow,dataCloseTS,dataHighTS,dataLowTS,dataLSA,dataLSB,dataConv,dataBl,dataLagSpan
},

raw=SemanticImport[FileNameJoin[{NotebookDirectory[],StringJoin[filenamE,".csv"]}]];
data=raw[starT;;enD,All];

{dataDates,dataClose,dataHigh,dataLow}=Normal@data[All,#]&/@{"Date","Close", "High","Low"};
{dataCloseTS,dataHighTS,dataLowTS}=TimeSeries[#,{dataDates}]&/@{dataClose,dataHigh,dataLow};

dataConv=TimeSeriesThread[Mean,{MovingMap[Max,dataHighTS,Quantity[clPeriod,"Events"]],MovingMap[Min,dataLowTS,Quantity[clPeriod,"Events"]]}];
dataBl=TimeSeriesThread[Mean,{MovingMap[Max,dataHighTS,Quantity[blPeriod,"Events"]],MovingMap[Min,dataLowTS,Quantity[blPeriod,"Events"]]}];
dataLSA=TimeSeriesShift[TimeSeriesThread[Mean,{TimeSeriesWindow[dataConv,{First[dataBl["Dates"]],Last[dataBl["Dates"]]}],dataBl}],{lagSpanPeriod,"BusinessDay"}];
dataLSB=TimeSeriesShift[TimeSeriesThread[Mean,{MovingMap[Max,dataHighTS,Quantity[leadSpanBPeriod,"Events"]],MovingMap[Min,dataLowTS,Quantity[leadSpanBPeriod,"Events"]]}],{lagSpanPeriod,"BusinessDay"}];
dataLagSpan=TimeSeriesShift[dataCloseTS,{{DatePlus[First[dataCloseTS["Dates"]],{-lagSpanPeriod,"BusinessDay"}]}}];

DateListPlot[
{dataLSA,dataLSB,dataConv,dataBl(*,dataLagSpan*),dataCloseTS},
PlotStyle->{Directive[Green,Thin],Directive[Red,Thin],Red,Blue(*Directive[Orange,Opacity[.7]],Directive[Gray,Opacity[.7]]*),Directive[Black,Opacity[.5]]},
Filling->{1->{{2},{LightRed,LightGreen}}},
PlotLegends->{"선행A","선행B","전환선","기준선"(*,"후행"*),"주가"},
(*GridLines->Automatic,*)PlotTheme->"Detailed",AspectRatio->1/1.8
]
]

고전적 방법: 이런 지표에서 적절한 전략을 만들어서 투자하는 방법

시계열 데이터의 구조

시계열 데이터는 다음과 같은 것들로 이뤄져 있다고 생각된다:

자기상관요인, 주기적 변동, 트렌드, 외인성, 화이트노이즈.

여기서 설명하는 내용은 대부분 Wolfram의 Tech Conference 등에서 발표된 내용에서 뽑은 것이다. (2012~2016) 이 내용을 설명한 분들은 다음과 같다: 
Gosia Konwerska, Arash Mahdian, Oleksandr Pavlyk, Anmol Bajracharya.

대표적 발표 내용은 2014년 Wolfram Technology Conference에서 행해진 Time Series Modeling Workshop을 들 수 있다.

예측을 위한 기본 가정: 정상성 Stationarity

확률과정이 정상적이란: 시간이 가도 그 통계적 특성이 변하지 않는다. 즉 확률분포가 변하지 않는다.

이것은 실제로 확인할 방법이 없다.

대신 다루기 쉬운 평균과 분산이 변하지 않음을 사용한다. (약 정상성 Weak stationarity)

실제로는 자기상관이 없음을 확인한다.

“ 定常”의 사전적 뜻: 일정하여 늘 한결같음.

백색소음의 지수이동평균 과정(WTC2014): (Subscript[y, t+1])==(1-α)(Subscript[y, t])+α (Subscript[ϵ, t+1]), (Subscript[y, 0])==(Subscript[ϵ, 0]).

dt=Table[{TemporalData[Drop[ExponentialMovingAverage[RandomVariate[NormalDistribution[],1100],α],100]],α},{α,{2.001,1.8,1,0.2,-0.001}}];

ListPlot[#[[1]],Filling->Axis,PlotLabel->If[0<#[[2]]<2,Row[{"stationary",", ",TraditionalForm[HoldForm[α]==#[[2]]]}],Row[{"non-stationary",", ",TraditionalForm[HoldForm[α]==#[[2]]]}]]]&/@dt

(*Table[ListPlot[Drop[ExponentialMovingAverage[RandomVariate[NormalDistribution[],1100],α],100],Filling->Axis,PlotLabel->If[0<α<2,Row[{"stationary",", ",TraditionalForm[HoldForm[α]==α]}],Row[{"non-stationary",", ",TraditionalForm[HoldForm[α]==α]}]]],{α,{2.001,1.8,1,0.2,-0.001}}]*)

공분산함수 (Subscript[γ, h ])==(n^(-1))(∑_(t=1)^(n-h))(Subscript[y, t])(Subscript[y, t-h]) (α==1인 경우는 순수한 백색소음이다).

자기상관계수 Autocorrelation

자기상관이 없으면 random walk에 가깝다고 생각된다.

여러 가지 경우들: 자기상관이 없다, 서서히 감소한다, 급격히 감소한다, 등등.

ListPlot[
CorrelationFunction[#[[1]],{0,20}],Filling->Axis,PlotRange->All,PlotLabel->TraditionalForm[HoldForm[α]==#[[2]]],PlotStyle->PointSize[Medium]]&/@dt

ListPlot[
CovarianceFunction[#[[1]],{0,20}],Filling->Axis,PlotRange->All,PlotLabel->TraditionalForm[HoldForm[α]==#[[2]]],PlotStyle->PointSize[Medium]]&/@dt

(*Table[ListPlot[CovarianceFunction[Drop[ExponentialMovingAverage[RandomVariate[NormalDistribution[],1100],α],100],{0,20}],Filling->Axis,PlotRange->All,PlotLabel->TraditionalForm[HoldForm[α]==α],PlotStyle->PointSize[Medium]],{α,{1.8,1.6,1,0.4,0.2}}]*)

계차수열을 잡으면

ListPlot[Differences[#[[1]]],Filling->Axis,PlotLabel->If[0<#[[2]]<2,Row[{"was stationary",", ",TraditionalForm[HoldForm[α]==#[[2]]]}],Row[{"was non-stationary",", ",TraditionalForm[HoldForm[α]==#[[2]]]}]]]&/@dt

dt2=Table[{TemporalData[Drop[ExponentialMovingAverage[RandomVariate[NormalDistribution[],1100],α],100]],α},{α,{1.8,1.6,1,0.4,0.2}}];

ListPlot[CorrelationFunction[#[[1]],{20}],Filling->Axis,PlotRange->All,PlotLabel->TraditionalForm[HoldForm[α]==#[[2]]],PlotStyle->PointSize[Medium]]&/@dt2

acf[data_,lmax_,clev_:0.95]:=Show[ListPlot[CorrelationFunction[data,{0,lmax}],Filling->Axis,PlotRange->{{0,lmax},All},PlotStyle->PointSize[Medium]],
Graphics[{Dashed,Line[{{0,#},{lmax,#}}]}]&/@((Quantile[NormalDistribution[],{((1-clev)/2),1-((1-clev)/2)}])/(Sqrt[data["PathLengths"][[1]]]))]

acf[#[[1]],20,.95]&/@dt2

주식 종가의 자기상관

ListLinePlot@sbux

ListPlot[CorrelationFunction[sbux["Values"],{40}],Filling->0]

acf[TemporalData@sbux["Values"],20,.95]

ListPlot[Differences@sbux,Filling->0,PlotStyle->PointSize[Medium]]

acf[Differences@TemporalData@sbux["Values"],20,.95]

편자기상관 Partial Auto-correlation

lag t와 t-k를 비교하는 상관계수를 계산할 때 그 사이의 값을 통해서 전달되는 상관성(일출溢出효과)을 제거하여, 둘 사이의 직접적인 관계만을 본다.

ListPlot[
PartialCorrelationFunction[#[[1]],{20}],Filling->Axis,PlotRange->All,PlotLabel->TraditionalForm[HoldForm[α]==#[[2]]],PlotStyle->PointSize[Medium]]&/@dt

ListPlot[
PartialCorrelationFunction[TemporalData@sbux["Values"],{20}],Filling->Axis,PlotRange->All,PlotLabel->"Partial Correlation of SBUX",PlotStyle->PointSize[Medium]]

계산의 예: 기계학습의 방법?

Ramesh Adhikari, “Selected Financial Applications”, The Mathematica Journal.

자기상관의 검증: AutocorrelationTest

timeSeriesData=TimeSeries[];(*WTI crude oil price*)

DateListPlot[timeSeriesData]

AutocorrelationTest[timeSeriesData,Automatic,"HypothesisTestData"]["TestDataTable",All]

Auto correlation이 없다.

UnitRootTest

AR 모델을 만족시키는 시계열이 전환함수의 분모에 unit root를 가지는지를 검증함.

unit root를 갖지 않으면 이 데이터에서 trend를 제거하면 정상시계열이 될 가능성이 있다.

UnitRootTest[timeSeriesData,Automatic,"HypothesisTestData"]["TestDataTable",All]

unit root를 갖는 AR 모델이라고 할 수 없다.

적절한 모형 찾기

TimeSeriesModelFit[timeSeriesData,Method->{Automatic,"SelectionCriterion"->#}]["CandidateSelectionTable"]&/@{"AIC","AICc","BIC","SBC"}

ARIMA 모델이 가장 적합하다고 하므로 이의 parameter를 추정한다.

timeseriesmodel=TimeSeriesModelFit[timeSeriesData,{"ARIMA",{1,1,0}}];

timeseriesmodel["ParameterTable"]

timeseriesmodel/@{"ACFPlot","PACFPlot","LjungBoxPlot"}

예측에 활용한다

Module[
{forecast,error,quantile,uband,lband},
forecast=TimeSeriesForecast[timeseriesmodel,{1,5}];
error=Sqrt@forecast@"MeanSquaredErrors";
quantile=Quantile[NormalDistribution[],1-1/2 (1-.95)];
uband=TimeSeriesThread[{1,quantile}.#&,{forecast,error}];
lband=TimeSeriesThread[{1,-quantile}.#&,{forecast,error}];
DateListPlot[
{
TimeSeriesWindow[timeSeriesData,{{2018,1,1},Last@timeSeriesData["Dates"]}],
forecast,
uband,
lband
},
PlotStyle->{Automatic,Automatic,Gray,Gray},GridLines->{Automatic,Automatic}, Filling->{3->{4}}]
]

계산의 예: Deep Neural Network

RNN의 기법을 쓴다: LSTM, GRU

통계적 의미를 파악하기는 어렵지만 유사한 기법을 사용한다.

Workflow의 예

데이터를 읽어들인다.

raw=FinancialData["SPY",{{2014,1,1},{2023,5,14}}];

(*Stationary하게 바꾼다*)
rawvalues=QuantityMagnitude@raw["Values"];
close=100Differences[Log[rawvalues]];
ln=Length@close;
y=close[[-ln+20-1;;]];
(*RNN의 input을 일정한 길이로 자른다.*)
input=Partition[close,20,1];(*Length/@{input,y}*)
(*train set과 test set으로 나눈다.*)
trainX=input[[;;ln-500]];trainy=y[[2;;ln-500+1]];
testX=input[[ln-500+1;;-2]];testy=y[[ln-500+2;;]];
(*Length/@{trainX,trainy,testX,testy}*)
trainXp=Transpose/@({#}&/@trainX);
testXp=Transpose/@({#}&/@testX);
Dimensions/@{trainXp,trainy,testXp,testy}

trainSet=MapThread[#1->#2&,{trainXp,#>0&/@trainy}];
testSet=MapThread[#1->#2&,{testXp,#>0&/@testy}];
Length/@{trainSet,testSet}

trainSet[[1]]

GRU 네트워크를 정의하고 Train한다.

Module[{lNet},
lNet=NetInitialize[#,All,RandomSeeding->1234]&@
NetChain[{
GatedRecurrentLayer[32,"Input"->{20,1}],
SequenceLastLayer[],
ElementwiseLayer[Ramp],
LinearLayer[16],
ElementwiseLayer[Ramp],
DropoutLayer[.5],
LinearLayer[],
ElementwiseLayer[LogisticSigmoid]},
"Output"->"Boolean"
];
trained=NetTrain[
lNet,
trainSet,
All,
LearningRate->Automatic,
ValidationSet->testSet,
Method->"ADAM",
BatchSize->32,
MaxTrainingRounds->100,
TargetDevice->"CPU"
]
];

결과를 확인한다.

trainedNet=trained["TrainedNet"];
NetMeasurements[trainedNet,testSet,"ConfusionMatrixPlot"]