3장 인자분석 - 주성분인자법(principal factor method) 사례

KNOU/3 다변량분석 2017. 3. 27. 19:59
반응형

### 5주차, 제3장 인자분석 ###

 

아래 자료는 무료 검진프로그램인 PHI를 개발하여 프로그램 유효서을 모니터링 한 자료이다. 11개 검진 항목, 128명의 자료가 이으며 R을 이용하여 주성분인자법(principal factor method)를 이용한 인자분석을 해 보자.

 

1. 자료 가져오기 및 요약 통계

> setwd('C:/Rwork')

> med.data = read.table('medFactor.txt',header=T)

> head(med.data)

lung muscle liver skeleton kidneys heart step stamina stretch blow urine

1   20     16    52       10      24    23   19      20      23   29    67

2   24     16    52        7      27    16   16      15      31   33    59

3   19     21    57       18      22    23   16      19      42   40    61

4   24     21    62       12      31    25   17      17      36   36    77

5   29     18    62       14      26    27   15      20      33   29    88

6   18     19    51       15      29    23   19      20      50   37    54

> boxplot(med.data)

 

2. 초기 인자분석 실행하기

> library(psych)

> library(GPArotation) #인자회전

> med.factor = principal(med.data,rotate='none') #rotate에 ‘none’, ‘varimax’,’quanrtmax’, ‘oblimin’ 등 방법을 사용

> names(med.factor)

[1] "values"       "rotation"     "n.obs"        "communality"  "loadings"     "fit"        

[7] "fit.off"      "fn"           "Call"         "uniquenesses" "complexity"   "chi"        

[13] "EPVAL"        "R2"           "objective"    "residual"     "rms"          "factors"    

[19] "dof"          "null.dof"     "null.model"   "criteria"     "STATISTIC"    "PVAL"       

[25] "weights"      "r.scores"     "Structure"    "scores"     

> med.factor$values #고유근

[1] 3.3791814 1.4827707 1.2506302 0.9804771 0.7688022 0.7330511 0.6403994 0.6221934 0.5283718 0.3519301

[11] 0.2621928

> plot(med.factor$values,type='b')

 

>> med.factor$values #고유근을 통해 세 번째 인자까지 고유근이 1 이상인 것을 알 수 있다. 스크리 그림에서는 4번째 인자다음부터 그래프의 기울기가 완만해지는 것을 볼 수 있다. 따라서 유효한 인자의 수는 3-4개로 생각할 수 있다.

 

3. 인자 회전

(1) 직교회전(orthogonal rotation)의 Varimax

> med.varimax = principal(med.data,nfactors=3,rotate='varimax')

> med.varimax

 

>> 공통성은 많은 문항(변수) 중에서 서로연관이 있는 일부 문항들을 선택하기 위하여 인자분석을 실시하는 경우 문항선택의 기준으로 이용된다. 예를 들면, 문항 100개를 이용하여 인자분석을 한 결과 공통성이 0.3 이하인 것이 40개라면 이 40개 문항은 다른 문항들과 공통점이 별로 없는 것으로 판단할 수 있다.


>> 고유분산은 u2 = 1 – h2(공통성) 공식으로 바로 구할 수 있다.lung의 경우 0.53 = 1 – 0.47


>> ss loading은 각 인자에 의해 설명되는 분산의 양을 나타낸다.

   RC1 2.39 = 0.662 + 0.112 + 0.782 + … + 0.252 + (-0.07)2


>> Proportion Var는 각 인자가 설명하는 총 분산의 비율을 의미한다. RC1은 총 분산의 22%, RC2가 19%, RC3는 14%, 아래 Cumulative var는 누적값. 세 인자에 의해 설명되어지는 변동은 총 변동의 56%.



>> 회전된 요인에 대한 변수들의 요인 적재값을 보면 첫 번째 인자는 lung, liver, kidneys, heart가 높은 값을 가지며, step은 상대적으로 높은 적재값을 가진다. (적재값이 뭘 의미하는지 잘 모르겠다 ㅠ)


두 번째 인자는 stamina, stretch, blow, urine의 값이 높다. 마지막으로 세번째 인자는 muscle과 skeleton에서 높은 값을 가진다. 따라서 첫 번째 인자는 생물의학(biomedical), 두번째 인자는 인체기능(performance), 세번째 인자는 근육골계통력(Muscular – skeletal strength)으로 특정 지을 수 있다.


>> 인자분석은 주성분분석과 달리 인자들의 결합으로 원변수들을 표현한다.

 Lung = 0.66 x factor1 + 0.12 x facotr2 + 0.16 x factor3 + error

 


인자점수 추정방법은 회귀분석, Barlett 등이 있는데, 아래는 회귀분석을 이용한 인자점수이다.

> head(med.varimax$scores)

             RC1        RC2         RC3

[1,] -0.11970907 -0.2574867 -0.74473196

[2,]  0.05696634 -0.4784611 -1.53199247

[3,] -0.59153602  0.8957933  1.52598533

[4,]  1.20919164  0.3179760 -0.42022213

[5,]  0.82291043  0.1513448 -0.02988744

[6,] -0.08606120  1.1451913  0.76290104 

 


인자분석의 시각화 몇가지 사례


> plot(med.varimax)

> plot(med.varimax$loadings[,c(1:2)])

> text(med.varimax$loadings[,c(1:2)],colnames(med.data))


> biplot(med.varimax$scores,med.varimax$loadings)

> fa.diagram(med.varimax)


> omega(med.data)

 





3. 인자 회전

(2) 사곽회전(Oblique rotation) - Oblimin

> med.oblimin = principal(med.data,nfactors=3,rotate='oblimin',scores=T,method='regression')

> med.oblimin

 


>> Oblimin을 이용한 인자분석에 대한 결과는 Varimax와 마찬가지로 세 인자에 대해 묶여지는 변수는 같으며 다만 인자 적재값이 차이가 나는 것을 볼 수 있다. Oblimin 방법을 이용한 결과에 대한 각 변수의 인자 모형은 아래와 같다.

Lung = 0.66 x factor1 + 0.01 x factor2 + 0.1 x facotr3

 

>> 회귀분석을 이용한 인자점수 함수식은 아래와 같다

Factor1 = 0.66 x lung + 0.02 x muscle + …+ (-0.09) x urine

 

인자점수는 아래와 같다.

> head(med.oblimin$scores)

TC1        TC2         TC3

[1,] -0.2325231 -0.2639103 -0.77351112

[2,] -0.1722964 -0.4641659 -1.54539480

[3,] -0.2852637  0.8308965  1.50251465

[4,]  1.1976234  0.4291411 -0.22225725

[5,]  0.8296036  0.2260879  0.09586509

[6,]  0.1766403  1.1289891  0.83993963 

 

 

4. 행렬도

> biplot(med.varimax)

 

 

참고문헌: 

[1] 다변량분석(김성수, 김현중, 정성석, 이용구 공저)

[2] R 프로그램에 기반한 다변량분석 및 데이터마이닝(이재길)


반응형

'KNOU > 3 다변량분석' 카테고리의 다른 글

2장 주성분분석 - heptathlon data / beer data  (0) 2017.03.27
Posted by 마르띤
,

2장 주성분분석 - heptathlon data / beer data

KNOU/3 다변량분석 2017. 3. 27. 19:13
반응형

[예제1] - heptathlon data



1. 자료 준비

> library(HSAUR2)

> library(tools)

> data("heptathlon")

> head(heptathlon)

> heptathlon은 R의 HSAUR2 패키지에 들어있는 자료로서, 1988년 서울 올림픽 여성 7종 경기에 대한 결과이다. Hurdles(110m 허들), highjump(높이뛰기), shot(포환던지기), run200m(200미터 달리기), longjump(멀리뛰기), javelin(창던지기), run800m(800미터 달리기),score(점수)로 구성되어 있다.


2. 자료 변형

> heptathlon$hurdles <- max(heptathlon$hurdles) - heptathlon$hurdles

> heptathlon$run200m <- max(heptathlon$run200m) - heptathlon$run200m

> heptathlon$run800m <- max(heptathlon$run800m) - heptathlon$run800m

>  작은 값일수록 좋은 점수이기 때문에 최대값에서 빼줌


 

3. 주성분분석

> library(stats)

> hep.data <- heptathlon[,-8] #주성분점수에서 예외

> heptathlon.pca <- princomp(hep.data,cor=T,scores=T) #cor = T, 상관계수 행령, F는 공분산 행렬, scores=T 주성분점수

> heptathlon.pca #변수가 7개이니, 주성분도 7개

Call:

  princomp(x = hep.data, cor = T, scores = T)

 

Standard deviations:

  Comp.1    Comp.2    Comp.3    Comp.4    Comp.5    Comp.6    Comp.7

2.1119364 1.0928497 0.7218131 0.6761411 0.4952441 0.2701029 0.2213617

 

7  variables and  25 observations.

> summary(heptathlon.pca)

> 주성분의 표준편차, 분산비율을 알수 있다. 첫 번째 주성분이 63.72%, 두 번째 주 성분이 17.06%의 분산 비율을 나타내고 있으며, 2개의 주성분이 총 변이의 80.8%의 정보를 가지고 있다.

 

> eig.val = heptathlon.pca$sdev^2 #주성분의 표준편차를 제곱

> eig.val #주성분 분산

Comp.1     Comp.2     Comp.3     Comp.4     Comp.5     Comp.6     Comp.7

4.46027516 1.19432056 0.52101413 0.45716683 0.24526674 0.07295558 0.04900101

> 제1주성분의 분산이 4.46, 제2주성분이 1.19, 제3주성분이 0.52로 유효한 주성분의 수는 2개임을 알 수 있다.

 

> par(mfrow=c(1,2))

> screeplot(heptathlon.pca)

> screeplot(heptathlon.pca,type='line',pch=19,main='Scree Plot of heptathlon.pca')

 

> 스크리 그림은 주성분의 고유값 크기 순으로 그린 것으로, 고유값이 1보다 큰 주성분값이 2개임을 알 수 있다.

 

> heptathlon.pca$loadings[,1:2]

             Comp.1      Comp.2

hurdles  -0.4528710  0.15792058

highjump -0.3771992  0.24807386

shot     -0.3630725 -0.28940743

run200m  -0.4078950 -0.26038545

longjump -0.4562318  0.05587394

javelin  -0.0754090 -0.84169212

run800m  -0.3749594  0.22448984

> 상관계수행렬을 이용한 경우 제1, 제2주성분은 다음과 같다.

PC1 = -0.453 x hurldes – 0.377 x highjump - ․․․ - 0.374 x rum800m

PC2 = 0.157 x hurldes + 0.248 x highjump + ․․․ + 0.224 x rum800m

 

여기서 제1주성분은 jvelin(창던지기)를 제외한 모든 변수의 절대값이 큰 값을 가지는 것으로 볼 때, 전반적인 체력 정도를 나타내는 성분이라고 할 수 있고, 제2성분은 jvelin의 계수가 다른 변수에 비해 상대적으로 절대값이 큰 것을 볼 때, 창던지기와 밀접한 관련이 있는 성분으로 파악할 수 있다.

 


4. 주성분 점수 및 행렬도 (biplot)

> head(heptathlon.pca$scores[,1:2])

                         Comp.1      Comp.2

Joyner-Kersee (USA) -4.206435 -1.26802363

John (GDR)          -2.941619 -0.53452561

Behmer (GDR)        -2.704271 -0.69275901

Sablovskaite (URS)  -1.371052 -0.70655862

Choubenkova (URS)   -1.387050 -1.78931718

Schulz (GDR)        -1.065372  0.08104469

> biplot(heptathlon.pca,cex=0.7,col=c('red','blue'),main='Biplot')

 

> 첫번째 주성분을 x축, 두 번째 주성분을 y축에 위치. Highjump, run800m, hurldes, longjump가 서로 

가까운 곳에 위치하고 벡터의 방향(화살표) 또한 비슷하다. Ru200m, shot이 가깝게 위치하고 있고, javelin은 다른 변수들과 다른 방향에 위치하고 있다.





[예제2] - beer data


99명의 소비자가 맥주 구매에 중요하게 생각하는 요인에 대해 점수를 준 엑셀 자료이다. 독립변수는 7개, 케이스 수는 99개인 자료이다. 점수는 0~100으로 되어 있다.

상관관계가 아닌 공분산을 활용해보자.

 

1. 자료 읽기

> beer<-read.csv('beer.csv',header=T)

> head(beer)

   cost size alcohol reputat color aroma taste

1   10   15      20      85    40    30    50

2  100   70      50      30    75    60    80

3   65   30      35      80    80    60    90

4    0    0      20      30    80    90   100

5   10   25      10     100    50    40    60

6   25   35      30      40    45    30    65

 

 

2. 상관계수행렬 및 산점도 행렬 보기

> round(cor(beer),2)

cost  size alcohol reputat color aroma taste

cost     1.00  0.88    0.88   -0.17  0.32 -0.03  0.05

size     0.88  1.00    0.82   -0.06  0.01 -0.29 -0.31

alcohol  0.88  0.82    1.00   -0.36  0.40  0.10  0.06

reputat -0.17 -0.06   -0.36    1.00 -0.52 -0.52 -0.63

color    0.32  0.01    0.40   -0.52  1.00  0.82  0.80

aroma   -0.03 -0.29    0.10   -0.52  0.82  1.00  0.87

taste    0.05 -0.31    0.06   -0.63  0.80  0.87  1.00

> 맥주의 가격(cost),크기(size), 알코올(alcohol)함량의 상관계수가 높아 서로 양의 상관관계가 있을 것으로 보인다. 또한 색(color), 향(aroma), taste(맛)이 서로 상관계수가 높게 나타났다. 반면 평판(reputation)은 비교적 다른 변수들과 상관계수가 높지 않으며, 다른 변수들과 모두 음의 상관 계수를 나타냈다.

 

> plot(beer[1:3])

 

 


3. 주성분분석 실행하기

> library(stats)

> beer.pca = princomp(beer,cor=F, scores=T) #cor = F 공분산

> beer.pca

Call:

  princomp(x = beer, cor = F, scores = T)

 

Standard deviations:

  Comp.1    Comp.2    Comp.3    Comp.4    Comp.5    Comp.6    Comp.7

39.117489 34.783073 17.974805  8.456192  6.479111  4.381410  2.579257

 

7  variables and  99 observations.

> summary(beer.pca)

 

 

> screeplot(beer.pca,type='l',pch=19)

> beer.pca$loadings[,1:3] #주성분계수

Comp.1      Comp.2     Comp.3

cost     0.7344197  0.31868378 -0.1902272

size     0.3942553  0.34017182  0.1274959

alcohol  0.2831873  0.06888179  0.0370581

reputat -0.3356901  0.49057105 -0.7861697

color    0.2657171 -0.34379889 -0.3862210

aroma    0.1398211 -0.48510439 -0.3693624

taste    0.1488354 -0.42871339 -0.2062207

> screeplot을 보면 4번째 주성분부터 완만해짐로 유효한 주성분을 3개로 판단하여, 3개의 주성분 계수를 보자. 제1주성분을 보면 cost, size, alcohol의 점수가 높다. 비용, 크기, 알코올이 서로 연관이 있고, color, aroma, taste의 주성분점수가 상대적으로 낮게 측정되어 있다. 색, 향, 맛이 서로 연관되어 있다.

PC1 = 0.734 x cost + 0.397 x size + … + 0.148 x taste

 


4. 주성분 분석 점수 및 행렬도

> head(beer.pca$scores[,1:3])

Comp.1    Comp.2     Comp.3

[1,] -41.435562  40.03360   4.340612

[2,]  91.264626  23.06214   7.798284

[3,]  31.574345  15.79053 -34.501270

[4,]  -9.770913 -59.53112  -0.351844

[5,] -39.816498  37.52890 -16.165598

[6,]  -1.035114  22.08072  34.760924

> biplot(beer.pca,col=c('red','blue'),main = 'Biplot of Beer.pca')

 

 


참고 문헌

[1] R을 이용한 다변량 분석, 김성수 김현중 정성석 이용구 공저

반응형

'KNOU > 3 다변량분석' 카테고리의 다른 글

3장 인자분석 - 주성분인자법(principal factor method) 사례  (0) 2017.03.27
Posted by 마르띤
,

제3장 예측데이터의 기초분석

KNOU/3 예측방법론 2017. 3. 22. 20:05
반응형

### 3,4주차, 제3장 예측데이터의 기초 분석 ###

 

3.1 시계열의 요약

> gdp <- read.csv("gdpq.csv", header = TRUE)

> head(gdp)

gdp   gdpsa

1 12807.5 14934.6

2 15732.5 15246.5

3 14621.6 15472.1

4 18689.5 16198.0

5 14512.5 16664.4

6 17709.4 17180.3

> gdp_o  <- ts(gdp[,1]/1000, start=1970, frequency=4)

> gdp_o

Qtr1     Qtr2     Qtr3     Qtr4

1970  12.8075  15.7325  14.6216  18.6895

1971  14.5125  17.7094  16.5814  19.5058

1972  15.2014  18.4940  17.5132  21.5480

1973  17.1063  20.9626  20.6994  24.7473

1974  19.7445  23.2128  22.2597  26.1339

1975  20.0790  24.6205  23.9055  29.4487

1976  22.5622  28.2495  27.2866  33.1492

1977  24.5630  30.6101  30.7895  38.4318

1978  28.1520  34.5775  33.4854  40.9892

1979  32.1360  38.2981  36.1002  42.1807

1980  32.0372  36.8905  36.8476  40.1279

1981  33.0530  38.8709  39.6106  45.1715

1982  36.2445  41.7473  43.0336  48.6735

1983  40.6951  46.9821  49.2677  53.4268

1984  46.3244  52.1396  53.7155  56.9611

1985  49.3720  55.4260  57.4681  62.4993

1986  54.3321  61.6800  65.9805  70.2836

1987  62.3064  70.6516  72.8342  77.4281

1988  72.3188  76.2441  80.4584  87.2238

1989  75.6243  81.9830  85.7357  94.2551

1990  82.9977  90.4118  94.2292 101.3471

1991  91.5677 100.2760 102.4124 110.5688

1992  98.9458 107.5298 107.0952 114.5932

1993 103.2117 113.8639 115.0267 123.1619

1994 112.7728 123.0717 124.0430 135.3119

1995 123.1052 134.6932 136.4655 145.1600

1996 132.1325 144.3954 146.0147 155.6439

1997 139.4080 154.6306 155.6320 161.8583

1998 134.4830 143.3159 144.5975 154.1904

1999 143.0859 159.4288 161.7758 174.1675

2000 161.1069 174.1233 176.7639 182.6340

2001 167.9311 181.3261 181.9274 191.0444

2002 179.0416 194.0302 194.2672 206.5295

2003 185.3180 197.5801 198.0618 214.5984

2004 194.9231 209.3292 207.5780 220.4750

2005 200.2166 216.3471 216.9831 231.6941

2006 212.4693 227.4156 227.9013 242.2628

2007 221.9311 239.4580 239.1251 256.0003

2008 234.1789 249.8902 246.9623 247.4675

2009 224.3595 244.7084 249.5196 263.0376

2010 243.8225 263.1901 260.7457 275.9080

2011 254.2480 272.3620 270.2561 285.2294

2012 261.4804 278.8836 274.4844 289.3663

2013 265.2912 285.3123 283.5950 300.6547

> gdp_sa <- ts(gdp[,2]/1000, start=1970, frequency=4)

> gdp_sa

Qtr1     Qtr2     Qtr3     Qtr4

1970  14.9346  15.2465  15.4721  16.1980

1971  16.6644  17.1803  17.3029  17.1614

1972  17.5141  17.9938  18.2651  18.9837

1973  19.7739  20.4964  21.4395  21.8057

1974  22.8617  22.6609  22.9105  22.9179

1975  23.5435  24.1886  24.7458  25.5757

1976  26.7061  27.7776  28.2670  28.4967

1977  28.8486  30.2254  31.9596  33.3608

1978  32.9248  33.9149  34.2955  36.0687

1979  37.5335  37.7262  36.7787  36.6767

1980  36.6760  36.2477  36.8647  36.1147

1981  37.2882  38.7201  39.7069  40.9907

1982  40.9807  41.7302  42.6719  44.3161

1983  45.5381  46.9978  48.6324  49.2035

1984  51.1110  52.0170  52.8702  53.1425

1985  54.4390  55.6003  56.4163  58.3099

1986  59.6636  62.0513  64.7921  65.7692

1987  67.5395  70.8755  71.8664  72.9389

1988  77.9135  76.8991  79.7696  81.6630

1989  81.4808  82.6901  85.2974  88.1298

1990  89.1127  90.7379  93.9483  95.1870

1991  98.3701 100.0831 102.2049 104.1669

1992 105.8722 106.9636 106.9913 108.3369

1993 110.1757 112.9698 115.0582 117.0605

1994 120.1321 122.1115 124.0878 128.8680

1995 131.3819 133.6788 136.2208 138.1425

1996 141.2176 143.2933 145.5552 148.1203

1997 149.0865 153.3882 154.8591 154.1950

1998 143.4587 141.9956 143.8028 147.3296

1999 151.7891 157.9815 161.6553 167.0321

2000 170.3475 172.5127 176.8317 174.9363

2001 177.4443 179.7859 182.2794 182.7193

2002 189.0423 192.4812 195.1160 197.2289

2003 195.8857 196.0955 199.1527 204.4243

2004 206.1241 207.7387 208.5458 209.8967

2005 211.6290 214.7745 218.1001 220.7373

2006 224.3165 225.6878 229.0608 230.9839

2007 234.2216 237.6344 240.2943 244.3642

2008 246.6202 247.5462 247.8380 236.4944

2009 236.6559 242.6501 250.8293 251.4898

2010 256.9193 260.5657 262.1134 264.0678

2011 267.5257 269.7017 271.9174 272.9508

2012 275.1821 276.0057 276.1202 276.9067

2013 279.2230 282.3449 285.3485 287.9793

> gdp_gr <- ts((gdp_sa-lag(gdp_sa,-1))/lag(gdp_sa,-1)*100,  start=c(1970,2), frequency=4)

> gdp_sa-lag(gdp_sa,-1)

Qtr1     Qtr2     Qtr3     Qtr4

1970            0.3119   0.2256   0.7259

1971   0.4664   0.5159   0.1226  -0.1415

1972   0.3527   0.4797   0.2713   0.7186

1973   0.7902   0.7225   0.9431   0.3662

1974   1.0560  -0.2008   0.2496   0.0074

1975   0.6256   0.6451   0.5572   0.8299

1976   1.1304   1.0715   0.4894   0.2297

1977   0.3519   1.3768   1.7342   1.4012

1978  -0.4360   0.9901   0.3806   1.7732

1979   1.4648   0.1927  -0.9475  -0.1020

1980  -0.0007  -0.4283   0.6170  -0.7500

1981   1.1735   1.4319   0.9868   1.2838

1982  -0.0100   0.7495   0.9417   1.6442

1983   1.2220   1.4597   1.6346   0.5711

1984   1.9075   0.9060   0.8532   0.2723

1985   1.2965   1.1613   0.8160   1.8936

1986   1.3537   2.3877   2.7408   0.9771

1987   1.7703   3.3360   0.9909   1.0725

1988   4.9746  -1.0144   2.8705   1.8934

1989  -0.1822   1.2093   2.6073   2.8324

1990   0.9829   1.6252   3.2104   1.2387

1991   3.1831   1.7130   2.1218   1.9620

1992   1.7053   1.0914   0.0277   1.3456

1993   1.8388   2.7941   2.0884   2.0023

1994   3.0716   1.9794   1.9763   4.7802

1995   2.5139   2.2969   2.5420   1.9217

1996   3.0751   2.0757   2.2619   2.5651

1997   0.9662   4.3017   1.4709  -0.6641

1998 -10.7363  -1.4631   1.8072   3.5268

1999   4.4595   6.1924   3.6738   5.3768

2000   3.3154   2.1652   4.3190  -1.8954

2001   2.5080   2.3416   2.4935   0.4399

2002   6.3230   3.4389   2.6348   2.1129

2003  -1.3432   0.2098   3.0572   5.2716

2004   1.6998   1.6146   0.8071   1.3509

2005   1.7323   3.1455   3.3256   2.6372

2006   3.5792   1.3713   3.3730   1.9231

2007   3.2377   3.4128   2.6599   4.0699

2008   2.2560   0.9260   0.2918 -11.3436

2009   0.1615   5.9942   8.1792   0.6605

2010   5.4295   3.6464   1.5477   1.9544

2011   3.4579   2.1760   2.2157   1.0334

2012   2.2313   0.8236   0.1145   0.7865

2013   2.3163   3.1219   3.0036   2.6308

> gdp_gr

Qtr1         Qtr2         Qtr3         Qtr4

1970               2.088438927  1.479683862  4.691670814

1971  2.879367823  3.095821032  0.713608028 -0.817781990

1972  2.055193632  2.738936057  1.507741555  3.934279035

1973  4.162518371  3.653806280  4.601295837  1.708062222

1974  4.842770468 -0.878324884  1.101456694  0.032299601

1975  2.729743999  2.740034404  2.303564489  3.353700426

1976  4.419820376  4.012191971  1.761851276  0.812608342

1977  1.234879828  4.772501959  5.737558477  4.384285160

1978 -1.306923095  3.007155700  1.122220617  5.170357627

1979  4.061138882  0.513408022 -2.511517195 -0.277334435

1980 -0.001908569 -1.167793653  1.702176966 -2.034466576

1981  3.249369370  3.840088822  2.548547137  3.233191209

1982 -0.024395778  1.828909706  2.256639077  3.853121141

1983  2.757462863  3.205447746  3.478035142  1.174320001

1984  3.876756735  1.772612549  1.640233001  0.515034935

1985  2.439666933  2.133213321  1.467617980  3.356476763

1986  2.321561176  4.001937530  4.416990458  1.508054223

1987  2.691685470  4.939331798  1.398085375  1.492352476

1988  6.820228986 -1.301956657  3.732813518  2.373585927

1989 -0.223112058  1.484153322  3.153098134  3.320617041

1990  1.115286770  1.823758005  3.538102601  1.318491128

1991  3.344049082  1.741382798  2.120038248  1.919673127

1992  1.637084333  1.030865515  0.025896660  1.257672353

1993  1.697297966  2.536040161  1.848635653  1.740249717

1994  2.623942320  1.647686172  1.618438886  3.852272343

1995  1.950755812  1.748262127  1.901573024  1.410724353

1996  2.226034711  1.469859281  1.578510649  1.762286748

1997  0.652307618  2.885371915  0.958939475 -0.428841444

1998 -6.962806836 -1.019875407  1.272715493  2.452525264

1999  3.026886654  4.079607824  2.325462159  3.326089525

2000  1.984887935  1.271048885  2.503583794 -1.071866639

2001  1.433664711  1.319625370  1.386927451  0.241332811

2002  3.460499247  1.819116674  1.368860959  1.082894278

2003 -0.681036096  0.107103275  1.559036286  2.647014075

2004  0.831505843  0.783314518  0.388516921  0.647771377

2005  0.825310736  1.486327488  1.548414733  1.209169551

2006  1.621474939  0.611323732  1.494542461  0.839558755

2007  1.401699426  1.457081670  1.119324475  1.693714749

2008  0.923212156  0.375476137  0.117876986 -4.577022087

2009  0.068289143  2.532875791  3.370779571  0.263326493

2010  2.158934478  1.419278349  0.593976874  0.745631471

2011  1.309474309  0.813379799  0.821537276  0.380041880

2012  0.817473332  0.299292723  0.041484650  0.284839718

2013  0.836491136  1.118066921  1.063805296  0.921960340

> mean(gdp_gr);median(gdp_gr);var(gdp_gr);sd(gdp_gr)

[1] 1.720042

[1] 1.618439

[1] 2.977614

[1] 1.725576

> (1+mean(gdp_gr)/100)^4-1 #연평균 7.05% 성장

[1] 0.07059725

> plot(gdp_gr, ylab="경제성장률(전 분기대비)", xlab="연도", col="steelblue", main="경제성장률 추이")

> abline(h=0, lty=2, col="gray")

3.2 시계열의 분포

> hist(gdp_gr, breaks = 12, col = "lightblue", border = "black", freq=FALSE, main="경제성장률의 확률분포", xlab="", xlim=c(-10, 10))

> lines(density(gdp_gr)) #꼬리부분에 큰 값을 가지고 있다. IMF



> shapiro.test(gdp_gr) #p-value = 4.285e-06, 경제성장률은 정규분포를 따른다는 귀무가설을 기각,

 

 

Shapiro-Wilk normality test

 

data:  gdp_gr

W = 0.94723, p-value = 4.285e-06

 

 

3.3 자기 상관

시계열은 시간의 전개에 따른 시간영역 정보와 시간에 따라 반복되는 주파수 영역 정보를 가지고 있다. 시간 영역 정보를 파악하려면 시계열의 과거, 현재, 미래를 연결하는 구조를 이해해야 한다. 시간에 따른 시계열의 연결 구조를 이해하는 데에는 표본 자기 상관계수, 표본 부분자기 상관계수 등이 이용된다.

> set.seed(1)

> nn=length(gdp_o)

> wn=ts(rnorm(nn), start=1970, frequency=4)

> wn

Qtr1         Qtr2         Qtr3         Qtr4

1970 -0.626453811  0.183643324 -0.835628612  1.595280802

1971  0.329507772 -0.820468384  0.487429052  0.738324705

1972  0.575781352 -0.305388387  1.511781168  0.389843236

1973 -0.621240581 -2.214699887  1.124930918 -0.044933609

1974 -0.016190263  0.943836211  0.821221195  0.593901321

1975  0.918977372  0.782136301  0.074564983 -1.989351696

1976  0.619825748 -0.056128740 -0.155795507 -1.470752384

1977 -0.478150055  0.417941560  1.358679552 -0.102787727

1978  0.387671612 -0.053805041 -1.377059557 -0.414994563

1979 -0.394289954 -0.059313397  1.100025372  0.763175748

1980 -0.164523596 -0.253361680  0.696963375  0.556663199

1981 -0.688755695 -0.707495157  0.364581962  0.768532925

1982 -0.112346212  0.881107726  0.398105880 -0.612026393

1983  0.341119691 -1.129363096  1.433023702  1.980399899

1984 -0.367221476 -1.044134626  0.569719627 -0.135054604

1985  2.401617761 -0.039240003  0.689739362  0.028002159

1986 -0.743273209  0.188792300 -1.804958629  1.465554862

1987  0.153253338  2.172611670  0.475509529 -0.709946431

1988  0.610726353 -0.934097632 -1.253633400  0.291446236

1989 -0.443291873  0.001105352  0.074341324 -0.589520946

1990 -0.568668733 -0.135178615  1.178086997 -1.523566800

1991  0.593946188  0.332950371  1.063099837 -0.304183924

1992  0.370018810  0.267098791 -0.542520031  1.207867806

1993  1.160402616  0.700213650  1.586833455  0.558486426

1994 -1.276592208 -0.573265414 -1.224612615 -0.473400636

1995 -0.620366677  0.042115873 -0.910921649  0.158028772

1996 -0.654584644  1.767287269  0.716707476  0.910174229

1997  0.384185358  1.682176081 -0.635736454 -0.461644730

1998  1.432282239 -0.650696353 -0.207380744 -0.392807929

1999 -0.319992869 -0.279113303  0.494188331 -0.177330482

2000 -0.505957462  1.343038825 -0.214579409 -0.179556530

2001 -0.100190741  0.712666307 -0.073564404 -0.037634171

2002 -0.681660479 -0.324270272  0.060160440 -0.588894486

2003  0.531496193 -1.518394082  0.306557861 -1.536449824

2004 -0.300976127 -0.528279904 -0.652094781 -0.056896778

2005 -1.914359426  1.176583312 -1.664972436 -0.463530401

2006 -1.115920105 -0.750819001  2.087166546  0.017395620

2007 -1.286300530 -1.640605534  0.450187101 -0.018559833

2008 -0.318068375 -0.929362147 -1.487460310 -1.075192297

2009  1.000028804 -0.621266695 -1.384426847  1.869290622

2010  0.425100377 -0.238647101  1.058483049  0.886422651

2011 -0.619243048  2.206102465 -0.255027030 -1.424494650

2012 -0.144399602  0.207538339  2.307978399  0.105802368

2013  0.456998805 -0.077152935 -0.334000842 -0.034726028

> par(mfrow=c(2,1))

> plot(wn, main="백색잡음계열", xlab="", ylab="", col="steelblue")

> abline(h=0, lty=2, col="gray")

> acf(wn, main="ACF 상관도표", col="steelblue")


 

>> 표준정규분포를 따르는 백색잡음계열의 상관도표는 위와 같다. 이를 보면 시차 0을 제외하고는 모두 점선 안의 작은 값을 갖는다. 시차의 1,2,3은 1,2,3년을 의미하며, 4개의 시차가 모여 1년이 된다. 백색잡음계열은 예측할 수 없는 계열이다. 시차 0을 제외한 모든 자기상관계수값이 점선 내에 있다. 따라서 백색잡음의 각 시차의 자기상관계수값이 5% 유의수준에서 0과 다르다고 할 수 없다.

 

 

 

> nn=length(gdp_sa)

> sin = ts(sin(1:nn/nn*12*pi), start=1970, frequency=4)

> par(mfrow=c(2,1))

> plot(sin, main="sin 계열", xlab="", ylab="", col="steelblue")

> abline(h=0, lty=2, col="gray")

> acf(sin, main="ACF 상관도표", col="steelblue")


>> sin 커브를 가지는 시계열의 상관도표를 구해보면 위와 같다. Sin 시계열의 자기상관계수는 sin커브와 같이 일정한 주기를 가지고 움직이는 것으로 나타났다.

 

> plot(gdp_o, main="GDP 원계열", xlab="", ylab="", col="steelblue")

> lines(gdp_sa,col='red')

> acf(gdp_o, main="ACF 상관도표", col="steelblue")

>> GDP원계열의 상관도표를 그려보면 위와 같다. 자기상관계수값이 서서히 감소하는 것으로 나타났다. 동 시계열은 강한 추세가 있는 불안정시계열임을 나타낸다. 강한 추세로 인해 GDP의 계절변동을 상관도표로 파악하기는 어렵다.

 

> plot(diff(log(gdp_o)), main="GDP로그차분", xlab="", ylab="", col="steelblue")

> acf(diff(log(gdp_o)), main="ACF 상관도표", col="steelblue")

>> GDP 원계열을 로그차분한 계열의 상관도표는 위와 같다. 이를 보면 자기상관계수값이 주기적으로 움직이는 것으로 나타났다. 이는 동 시계열의 로그차분을 통해 추세를 제거했지만, 계절변동은 그대로 남아 있음을 의미한다.

 

> plot(diff(log(gdp_o),4), main="GDP로그4차 차분", xlab="", ylab="", col="steelblue")

> acf(diff(log(gdp_o),4), main="ACF 상관도표 ", col="steelblue")

>>GDP 원계열을 로그 4차 차분한 계열의 상관도표를 그려보면 위와 같다. 자기상관계수값이 3차 이후 급격히 하락하는 것으로 나타났다. 이는 동 시계열의 로그 4차 차분을 통해 추세 및 계절 변동이 제거되었음을 의미한다.

 

 

륭(Ljung) – 박스(Box) 검정을 통해, 시계열에 m차까지의 자기상관관계가 존재하지 않는다는 귀무가설을 검정하는 것이다.

 - 귀무가설: H0 : p(1) = p(2) = … = p(m) = 0 잔차는 랜덤

> Box.test(wn,lag=8, type="Ljung")

 

Box-Ljung test

 

data:  wn

X-squared = 3.16, df = 8, p-value = 0.9239

 

> Box.test(sin,lag=8, type="Ljung")

 

Box-Ljung test

 

data:  sin

X-squared = 579.9, df = 8, p-value < 2.2e-16

 

> Box.test(gdp_o,lag=8, type="Ljung")

 

Box-Ljung test

 

data:  gdp_o

X-squared = 1249.8, df = 8, p-value < 2.2e-16

 

> Box.test(diff(log(gdp_o)),lag=8, type="Ljung")

 

Box-Ljung test

 

data:  diff(log(gdp_o))

X-squared = 756.05, df = 8, p-value < 2.2e-16

 

> Box.test(diff(log(gdp_o),4),lag=8, type="Ljung")

 

Box-Ljung test

 

data:  diff(log(gdp_o), 4)

X-squared = 205.63, df = 8, p-value < 2.2e-16

 

계열명

X2 검정통계량 값

유의확률

백색잡음계열

3.16

0.9239

sin계열

579.9

< 2.2e-16

GDP원계열

1249.8

< 2.2e-16

GDP 로그차분계열

756.05

< 2.2e-16

GDP 로그 4차 차분계열

205.63

< 2.2e-16

>> 백색잡음계열의 경우 p-value가 0.05보다 크므로 시계열에 8차까지의 자기상관관계가 존재하지 않는다는 귀무가설을 기각하지 못하지만, 나머지 계열은 모두 작은 값이어서 8차까지의 자기상관관계가 존재하지 않는다는 귀무가설을 기각하낟.

 

 

3.4 부분 자기 상관

표본 부분 자기 상관계수(sample partial autocorrelation coefficietns). 앞뒤 시차 영향을 제외한 순수한 영역

> gdp <- read.csv("gdpq.csv", header = TRUE)

> gdp_o  <- ts(gdp[,1]/1000, start=1970, frequency=4)

> gdp_sa <- ts(gdp[,2]/1000, start=1970, frequency=4)

> gdp_gr <- ts((gdp_sa-lag(gdp_sa,-1))/lag(gdp_sa,-1)*100,  start=c(1970,2), frequency=4)

> par(mfrow=c(2,1))

> plot(wn, main="", xlab="백색잡음계열", ylab="", col="steelblue")

> abline(h=0, lty=2, col="gray")

> pacf(wn, main="Partial ACF 부분상관도표", col="steelblue") #main="", xlab="", ylab="", col="steelblue")

> abline(h=0, lty=2, col="gray")


>> 모두 점선 안(1.96 / root N)의 작은 값을 가지고 있다.

 

> plot(sin, main="sin 계열", xlab="", ylab="", col="steelblue")

> abline(h=0, lty=2, col="gray")

> pacf(sin, main="Partial ACF 부분상관도표", col="steelblue")


>> sin 계열의 부분자기상관계수는 앞의 2기에서 유의하게 큰 값을 가지고 있다.

 

> plot(gdp_o, main="GDP 원계열", xlab="", ylab="", col="steelblue")

> lines(gdp_sa,col='red')

> pacf(gdp_o, main="Partial ACF 부분상관도표", col="steelblue")


>> 부분자기상관계수값은 1차에서 1에 가까운 큰 값을 가지고 있다. 1차의 큰 값은 추세가 있음을, 주기적으로 크게 나타난 부분자기상관계수는 계절변동이 있음을 의미한다.

 

 

> plot(diff(log(gdp_o)), main="GDP 로그차분", xlab="", ylab="", col="steelblue")

> pacf(diff(log(gdp_o)), main="Partial ACF 부분상관도표 ", col="steelblue")


>> 부분자기상관계수값이 1차의 값은 작고 주기적으로 움직이는 것으로 나타났다. 이는 동 시계열의 로그차분을 통해 추세를 제거했지만, 계절 변동은 그대로 남아있음을 의미한다.



 

> plot(diff(log(gdp_o),4), main="GDP 로그 4차 차분", xlab="", ylab="", col="steelblue")

> abline(h=0,lty=2,col='grey')

> pacf(diff(log(gdp_o),4), main="Partial ACF 부분상관도표 ", col="steelblue")

>> GDP 원계열을 로그 4차차분한 계열의 부분상관도표를 그려보면 그림과 같다. 이를 살펴보면 부분자기상관계수값은 3차 이후 작은값을 가지는 것으로 나타났다. 이는 동시계열의 로그 4차 차분을 통해 추세 및 계절변동이 제거되었음을 의미한다.

 

> sin12 = ts(sin(1:nn/nn*12*pi), start=1970, frequency=4)

> sin36 = ts(sin(1:nn/nn*36*pi), start=1970, frequency=4)

> plot(cbind(sin12, sin36), main="sin12와 sin36 추이", xlab="", ylab="", col="steelblue")

> spectrum(cbind(sin12, sin36), spans=c(3,3), main="sin12와 sin36 함수 스펙트럼", col="steelblue")

>> 위의 그래프를 보면 첫번째 sin함수는 주기가 상대적으로 긴 함수, 즉 저주파 함수이고, 두번째 sin 함수는 상대적으로 주기가 짧은 함수, 즉 고주파 함수이다. 이를 평활화된 주기도(스펙트럼)을 구해보면 아래 그래프와 같은데 실선이 sin12, 점선이 sin36이다. 해당주파수(주기의 역수)의 값이 크게 나타난다. 저주파 sin12자료의 스펙트럼은 낮은 주파수에서 최고점이 나타나나, 고주파 sin36 자료의 스펙트럼은 높은 주파수에서 최고점이 나타난다.

 

> plot(sin12+sin36, main="sin12+sin36", xlab="", ylab="", col="steelblue")

> spectrum(sin12+sin36, spans=c(3,3), main="sin 12 + sin36 합성 스펙트럼", col="steelblue")

 

>> 두 sin 계열을 합성한 결과 스펙트럼은 다음과 같다. 두 계열 합성계열의 스펙트럼은 두 sin 계열에 해당되는 주파수에서 큰 값을 가지고 있어서 동 계열에 두 주파수에서의 변동요인이 있음을 알 수 있다.

 

> par(mfrow=c(2,1))

> plot(wn, main="백색잡음계열", xlab="", ylab="", col="steelblue")

> abline(h=0, lty=2, col="gray")

> wn=ts(rnorm(nn), start=1970, frequency=4)

> spectrum(wn, spans=c(3,3), main="스펙트럼", col=1:2)

>> 표준정규분포로부터 154개의 난수를 생성하여 보면 아무 특징이 없는 백색잡음계열이 생성된다. 백색잡음계열의 평활화한 스펙트럼은 그림과 같이 아무 특징이 없는데 y축 값을 보면 매우 작은값을 가지고 수평적으로 움직이고 있다. 백색잡음계열의 평활화한 스펙트럼은 이론적으로는 수평선이다.  다양한 peak가 나타난다는 것은 잔차가 random, 주파수에 특정 주기가 없다는 뜻이다.

 

> par(mfrow=c(2,1))

> plot(gdp_o, main="GDP 원계열", xlab="", ylab="")

> lines(gdp_sa, col=2)

> spectrum(cbind(gdp_o, gdp_sa), spans=c(3,3), main="스펙트럼(원계열+계정조정계열)", col=1:2)

 

>> GDP 원계열과 계절조정계열의 시계열도표와 스펙트럼은 다음과 같다. 원계열의 스펙트럼(실선)을 보면 저주파수와 계절주파수에서 큰 값을 가지는 것으로 나타났다. 여기서 저주파수는 추세변동요인의 존재를, 계절주파수의 큰 값은 계절변동요인의 존재를 의미한다. 반면 GDP 계절조정계열의 스펙트럼(점선)을 구해보면 큰 값이 사라짐을 알 수 있다. 즉 계절변동요인이 제거되었음을 말한다.

 

> dlgdp1 = diff(log(gdp_o))

> dlgdp2 = diff(log(gdp_o),4)

> dlgdp  = cbind(dlgdp1,dlgdp2)

> plot(dlgdp1, main="로그1차차분과 4차차분계", xlab="", ylab="", col="steelblue")

> lines(dlgdp2, col=2)

> spectrum( na.omit(cbind(dlgdp1,dlgdp2)), spans=c(3,3), main="스펙트럼(로그1차차분과 4차차분계열)", col=c("steelblue", "red"))

 

>> 로그1차 차분한 GDP 원계열(보라색)의 스펙트럼을 구해보면 아래 그림과 같다(보라색). 이를 보면 원계열에서의 스텍트럼에서 계절주파수의 큰 값은 남고 저주파수의 추세변동요인과 관련된 큰 값은 사라짐을 알 수 있다. 이는 차분이 시계열에서 추세변동을 제거하는 역할을 했음을 의미한다. 한편 로그 4차 차분 계열의 스펙트럼(점선)을 구해보면 추세변동요인과 계졀변동 요인에 해당되는 주파스의 밀도함수값이 작게 나타난반면, 순환변동에 해당하는 값이 상대적으로 크게 나타나 로그 4차 차분계열에서 순환변동을 파악할 수 있을 것으로 보인다.


출처: 예측방법론(이긍희, 이한식, 장영재 공저)

반응형

'KNOU > 3 예측방법론' 카테고리의 다른 글

제2장 예측데이터 - 시계열  (1) 2017.03.17
Posted by 마르띤
,

제2장 예측데이터 - 시계열

KNOU/3 예측방법론 2017. 3. 17. 19:04
반응형

### 2주차, 제2장 예측데이터 - 시계열 ###

 

> setwd("c:/Rwork/")

> gdp<-read.csv('gdpq.csv',header=T)

> head(gdp)

gdp   gdpsa #계절 조정

1 12807.5 14934.6

2 15732.5 15246.5

3 14621.6 15472.1

4 18689.5 16198.0

5 14512.5 16664.4

6 17709.4 17180.3

> gdp <- ts(read.csv("gdpq.csv", header = TRUE), start=1970, frequency=4)

> gdp

gdp    gdpsa

1970 Q1  12807.5  14934.6

1970 Q2  15732.5  15246.5

1970 Q3  14621.6  15472.1

1970 Q4  18689.5  16198.0

(…중간 생략…)

2013 Q1 265291.2 279223.0

2013 Q2 285312.3 282344.9

2013 Q3 283595.0 285348.5

2013 Q4 300654.7 287979.3

> plot(gdp[,1]/1000, ylab="GDP(조원)", xlab="연도", col="steelblue")

> lines(gdp[,2]/1000, col="red")


 

> ww=c(.5,1,1,1,.5) #5분기 이동평균

> gdpm5 = filter(gdp, sides=2, ww/sum(ww))

> dlgdp1 = diff(log(gdp)) #1차 차분

> dlgdp4 = diff(log(gdp),4) #4차차분

> plot(gdp, main="GDP", ylab=" ", xlab="YEAR") #gdp 원계열


> plot(gdpm5, main="GDPM5", ylab=" ", xlab="YEAR") #gdp m5  이동평균


> plot(dlgdp1, main="diff(log(GDP))", ylab=" ", xlab= "YEAR") #1차 차분


> plot(dlgdp4, main="diff(log(GDP),4)", ylab=" ", xlab= "YEAR") #4차 차분

 

 

 

 출처: 예측방법론(이긍희, 이한식, 장영재 공저, Knou press)

반응형

'KNOU > 3 예측방법론' 카테고리의 다른 글

제3장 예측데이터의 기초분석  (0) 2017.03.22
Posted by 마르띤
,

티스토리툴바