R version 2.12.0 (2010-10-15)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(96
+ ,6.08
+ ,54.7
+ ,1914
+ ,1005
+ ,2
+ ,89
+ ,5.73
+ ,54.2
+ ,1684
+ ,963
+ ,2
+ ,87
+ ,6.22
+ ,53
+ ,1902
+ ,1035
+ ,2
+ ,87
+ ,5.8
+ ,52.9
+ ,1860
+ ,1027
+ ,2
+ ,101
+ ,7.99
+ ,57.8
+ ,2264
+ ,1281
+ ,2
+ ,103
+ ,8.42
+ ,56.9
+ ,2216
+ ,1272
+ ,2
+ ,103
+ ,7.44
+ ,56.6
+ ,1866
+ ,1051
+ ,2
+ ,96
+ ,6.84
+ ,55.3
+ ,1850
+ ,1079
+ ,2
+ ,127
+ ,6.48
+ ,53.1
+ ,1743
+ ,1034
+ ,2
+ ,126
+ ,6.43
+ ,54.8
+ ,1709
+ ,1070
+ ,2
+ ,101
+ ,7.99
+ ,57.2
+ ,1689
+ ,1173
+ ,1
+ ,96
+ ,8.76
+ ,57.2
+ ,1806
+ ,1079
+ ,1
+ ,93
+ ,6.32
+ ,57.2
+ ,2136
+ ,1067
+ ,1
+ ,88
+ ,6.32
+ ,57.2
+ ,2018
+ ,1104
+ ,1
+ ,94
+ ,7.6
+ ,55.8
+ ,1966
+ ,1347
+ ,1
+ ,85
+ ,7.62
+ ,57.2
+ ,2154
+ ,1439
+ ,1
+ ,97
+ ,6.03
+ ,57.2
+ ,1767
+ ,1029
+ ,1
+ ,114
+ ,6.59
+ ,56.5
+ ,1827
+ ,1100
+ ,1
+ ,113
+ ,7.52
+ ,59.2
+ ,1773
+ ,1204
+ ,1
+ ,124
+ ,7.67
+ ,58.5
+ ,1971
+ ,1160
+ ,1
+ ,129
+ ,7.57
+ ,57.3
+ ,2067
+ ,1401
+ ,1
+ ,110
+ ,6.45
+ ,53.7
+ ,2221
+ ,1142
+ ,1
+ ,102
+ ,7.99
+ ,56.6
+ ,2151
+ ,1288
+ ,1
+ ,134
+ ,8.43
+ ,57.5
+ ,2018
+ ,979
+ ,1
+ ,119
+ ,7.02
+ ,55.5
+ ,1677
+ ,1104
+ ,2
+ ,139
+ ,5.21
+ ,55.7
+ ,2126
+ ,956
+ ,2
+ ,75
+ ,6.21
+ ,53.1
+ ,1841
+ ,1153
+ ,1
+ ,138
+ ,5.39
+ ,55.9
+ ,2062
+ ,1001
+ ,2
+ ,132
+ ,5.59
+ ,57.8
+ ,1809
+ ,1230
+ ,1
+ ,122
+ ,7.72
+ ,59
+ ,2326
+ ,1014
+ ,2
+ ,102
+ ,6.69
+ ,58.4
+ ,1871
+ ,1287
+ ,1
+ ,78
+ ,5.96
+ ,55.4
+ ,2108
+ ,1198
+ ,2
+ ,119
+ ,8.49
+ ,59.5
+ ,1753
+ ,1125
+ ,2
+ ,136
+ ,6.64
+ ,53
+ ,1942
+ ,1142
+ ,1
+ ,109
+ ,5.23
+ ,54.6
+ ,1661
+ ,1379
+ ,2
+ ,85
+ ,6.2
+ ,58.4
+ ,1782
+ ,1148
+ ,2
+ ,119
+ ,7.36
+ ,58.2
+ ,2206
+ ,1318
+ ,2
+ ,136
+ ,6.67
+ ,53.2
+ ,1710
+ ,1041
+ ,2
+ ,72
+ ,6.36
+ ,54.2
+ ,2114
+ ,1253
+ ,2
+ ,125
+ ,7.43
+ ,53.8
+ ,2339
+ ,1264
+ ,1
+ ,87
+ ,8.41
+ ,53.8
+ ,2310
+ ,953
+ ,1
+ ,106
+ ,7.15
+ ,57.3
+ ,2012
+ ,1049
+ ,1
+ ,99
+ ,5.36
+ ,53
+ ,2028
+ ,1392
+ ,2
+ ,123
+ ,7.39
+ ,52.1
+ ,1728
+ ,1135
+ ,1
+ ,99
+ ,5.63
+ ,52.7
+ ,2107
+ ,1450
+ ,1
+ ,88
+ ,8.47
+ ,55.5
+ ,1965
+ ,958
+ ,1
+ ,97
+ ,7.75
+ ,57.8
+ ,2113
+ ,1209
+ ,2
+ ,119
+ ,8.33
+ ,55.4
+ ,1917
+ ,1441
+ ,2
+ ,77
+ ,6
+ ,57.9
+ ,2036
+ ,994
+ ,1
+ ,128
+ ,5.45
+ ,55.2
+ ,1869
+ ,1149
+ ,1
+ ,100
+ ,8.28
+ ,58.5
+ ,1858
+ ,1204
+ ,1
+ ,116
+ ,5.6
+ ,53.4
+ ,2145
+ ,1414
+ ,2
+ ,76
+ ,7.38
+ ,58.6
+ ,2009
+ ,1339
+ ,2
+ ,76
+ ,7.99
+ ,53.5
+ ,2249
+ ,1255
+ ,1
+ ,100
+ ,6.83
+ ,53.3
+ ,1949
+ ,1189
+ ,2
+ ,105
+ ,5.64
+ ,53.4
+ ,2058
+ ,1298
+ ,2
+ ,120
+ ,8.43
+ ,57.2
+ ,1953
+ ,1167
+ ,2
+ ,97
+ ,7.38
+ ,54.2
+ ,2089
+ ,1290
+ ,2
+ ,95
+ ,6.55
+ ,55.7
+ ,2254
+ ,1057
+ ,2
+ ,101
+ ,5.71
+ ,59.2
+ ,1973
+ ,1018
+ ,1)
+ ,dim=c(6
+ ,60)
+ ,dimnames=list(c('IQ'
+ ,'CCMIDSA'
+ ,'HC'
+ ,'TOTSA'
+ ,'TOTVOL'
+ ,'SEX
')
+ ,1:60))
> y <- array(NA,dim=c(6,60),dimnames=list(c('IQ','CCMIDSA','HC','TOTSA','TOTVOL','SEX
'),1:60))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> library(lattice)
> library(lmtest)
Loading required package: zoo
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
IQ CCMIDSA HC TOTSA TOTVOL SEX\r
1 96 6.08 54.7 1914 1005 2
2 89 5.73 54.2 1684 963 2
3 87 6.22 53.0 1902 1035 2
4 87 5.80 52.9 1860 1027 2
5 101 7.99 57.8 2264 1281 2
6 103 8.42 56.9 2216 1272 2
7 103 7.44 56.6 1866 1051 2
8 96 6.84 55.3 1850 1079 2
9 127 6.48 53.1 1743 1034 2
10 126 6.43 54.8 1709 1070 2
11 101 7.99 57.2 1689 1173 1
12 96 8.76 57.2 1806 1079 1
13 93 6.32 57.2 2136 1067 1
14 88 6.32 57.2 2018 1104 1
15 94 7.60 55.8 1966 1347 1
16 85 7.62 57.2 2154 1439 1
17 97 6.03 57.2 1767 1029 1
18 114 6.59 56.5 1827 1100 1
19 113 7.52 59.2 1773 1204 1
20 124 7.67 58.5 1971 1160 1
21 129 7.57 57.3 2067 1401 1
22 110 6.45 53.7 2221 1142 1
23 102 7.99 56.6 2151 1288 1
24 134 8.43 57.5 2018 979 1
25 119 7.02 55.5 1677 1104 2
26 139 5.21 55.7 2126 956 2
27 75 6.21 53.1 1841 1153 1
28 138 5.39 55.9 2062 1001 2
29 132 5.59 57.8 1809 1230 1
30 122 7.72 59.0 2326 1014 2
31 102 6.69 58.4 1871 1287 1
32 78 5.96 55.4 2108 1198 2
33 119 8.49 59.5 1753 1125 2
34 136 6.64 53.0 1942 1142 1
35 109 5.23 54.6 1661 1379 2
36 85 6.20 58.4 1782 1148 2
37 119 7.36 58.2 2206 1318 2
38 136 6.67 53.2 1710 1041 2
39 72 6.36 54.2 2114 1253 2
40 125 7.43 53.8 2339 1264 1
41 87 8.41 53.8 2310 953 1
42 106 7.15 57.3 2012 1049 1
43 99 5.36 53.0 2028 1392 2
44 123 7.39 52.1 1728 1135 1
45 99 5.63 52.7 2107 1450 1
46 88 8.47 55.5 1965 958 1
47 97 7.75 57.8 2113 1209 2
48 119 8.33 55.4 1917 1441 2
49 77 6.00 57.9 2036 994 1
50 128 5.45 55.2 1869 1149 1
51 100 8.28 58.5 1858 1204 1
52 116 5.60 53.4 2145 1414 2
53 76 7.38 58.6 2009 1339 2
54 76 7.99 53.5 2249 1255 1
55 100 6.83 53.3 1949 1189 2
56 105 5.64 53.4 2058 1298 2
57 120 8.43 57.2 1953 1167 2
58 97 7.38 54.2 2089 1290 2
59 95 6.55 55.7 2254 1057 2
60 101 5.71 59.2 1973 1018 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) CCMIDSA HC TOTSA TOTVOL `SEX\r`
118.482164 0.181843 0.288059 -0.013640 -0.004619 1.099468
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-30.826 -11.735 -2.589 15.002 34.742
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 118.482164 75.511825 1.569 0.122
CCMIDSA 0.181843 2.619733 0.069 0.945
HC 0.288059 1.267239 0.227 0.821
TOTSA -0.013640 0.013631 -1.001 0.321
TOTVOL -0.004619 0.017720 -0.261 0.795
`SEX\r` 1.099468 4.935794 0.223 0.825
Residual standard error: 18.64 on 54 degrees of freedom
Multiple R-squared: 0.02445, Adjusted R-squared: -0.06588
F-statistic: 0.2706 on 5 and 54 DF, p-value: 0.9272
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.37747546 0.75495092 0.62252454
[2,] 0.27538757 0.55077514 0.72461243
[3,] 0.15724098 0.31448196 0.84275902
[4,] 0.09829701 0.19659401 0.90170299
[5,] 0.19932004 0.39864008 0.80067996
[6,] 0.13152092 0.26304184 0.86847908
[7,] 0.08830477 0.17660954 0.91169523
[8,] 0.06079622 0.12159243 0.93920378
[9,] 0.03666185 0.07332370 0.96333815
[10,] 0.04161304 0.08322608 0.95838696
[11,] 0.02383475 0.04766951 0.97616525
[12,] 0.03988366 0.07976731 0.96011634
[13,] 0.08646025 0.17292049 0.91353975
[14,] 0.10171407 0.20342814 0.89828593
[15,] 0.06734915 0.13469829 0.93265085
[16,] 0.10652381 0.21304762 0.89347619
[17,] 0.08323893 0.16647785 0.91676107
[18,] 0.20144123 0.40288246 0.79855877
[19,] 0.25383776 0.50767551 0.74616224
[20,] 0.35411149 0.70822298 0.64588851
[21,] 0.40276307 0.80552613 0.59723693
[22,] 0.47767801 0.95535602 0.52232199
[23,] 0.41019302 0.82038605 0.58980698
[24,] 0.46812876 0.93625753 0.53187124
[25,] 0.41274224 0.82548448 0.58725776
[26,] 0.60249012 0.79501976 0.39750988
[27,] 0.54155853 0.91688293 0.45844147
[28,] 0.59188963 0.81622075 0.40811037
[29,] 0.65813164 0.68373672 0.34186836
[30,] 0.69884424 0.60231152 0.30115576
[31,] 0.79885905 0.40228190 0.20114095
[32,] 0.94524304 0.10951392 0.05475696
[33,] 0.93386394 0.13227212 0.06613606
[34,] 0.91766543 0.16466914 0.08233457
[35,] 0.90712631 0.18574738 0.09287369
[36,] 0.86391627 0.27216746 0.13608373
[37,] 0.79391921 0.41216159 0.20608079
[38,] 0.73604574 0.52790851 0.26395426
[39,] 0.65263025 0.69473950 0.34736975
[40,] 0.57301045 0.85397910 0.42698955
[41,] 0.62487790 0.75024420 0.37512210
[42,] 0.54058890 0.91882221 0.45941110
[43,] 0.39008527 0.78017053 0.60991473
> postscript(file="/var/www/rcomp/tmp/1nof71321784479.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/2k6sf1321784479.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/3x9pv1321784479.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/4lab51321784479.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/5vl651321784479.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 60
Frequency = 1
1 2 3 4 5 6
-10.79328689 -20.91692709 -19.35415005 -19.85882366 -0.98449753 0.50024782
7 8 9 10 11 12
-5.03015710 -11.63548077 18.39631464 16.61822724 -8.05433859 -12.03264222
13 14 15 16 17 18
-10.14303087 -16.58168833 -9.99797019 -16.41550723 -11.29915586 6.94705283
19 20 21 22 23 24
4.74400625 18.41592902 26.20252558 9.34742375 -0.04839387 28.37078977
25 26 27 28 29 30
9.02986177 34.74228468 -30.56865571 32.98682251 25.10940370 19.33127563
31 32 33 34 35 36
-4.15444901 -25.43532962 8.74399715 31.70883034 0.66668186 -24.02090088
37 38 39 40 41 42
16.39460905 26.91515923 -30.82649004 26.31354313 -13.69684737 0.90267055
43 44 45 46 47 48
-3.82996813 15.88031058 -1.34766344 -17.88031560 -7.33315454 13.65087545
49 50 51 52 53 54
-27.98773831 22.32807601 -7.03311489 14.70872358 -29.31441343 -23.97108520
55 56 57 58 59 60
-5.19901383 1.97888868 13.33954454 -6.18206542 -7.28885486 -5.05796482
> postscript(file="/var/www/rcomp/tmp/6gpyv1321784479.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -10.79328689 NA
1 -20.91692709 -10.79328689
2 -19.35415005 -20.91692709
3 -19.85882366 -19.35415005
4 -0.98449753 -19.85882366
5 0.50024782 -0.98449753
6 -5.03015710 0.50024782
7 -11.63548077 -5.03015710
8 18.39631464 -11.63548077
9 16.61822724 18.39631464
10 -8.05433859 16.61822724
11 -12.03264222 -8.05433859
12 -10.14303087 -12.03264222
13 -16.58168833 -10.14303087
14 -9.99797019 -16.58168833
15 -16.41550723 -9.99797019
16 -11.29915586 -16.41550723
17 6.94705283 -11.29915586
18 4.74400625 6.94705283
19 18.41592902 4.74400625
20 26.20252558 18.41592902
21 9.34742375 26.20252558
22 -0.04839387 9.34742375
23 28.37078977 -0.04839387
24 9.02986177 28.37078977
25 34.74228468 9.02986177
26 -30.56865571 34.74228468
27 32.98682251 -30.56865571
28 25.10940370 32.98682251
29 19.33127563 25.10940370
30 -4.15444901 19.33127563
31 -25.43532962 -4.15444901
32 8.74399715 -25.43532962
33 31.70883034 8.74399715
34 0.66668186 31.70883034
35 -24.02090088 0.66668186
36 16.39460905 -24.02090088
37 26.91515923 16.39460905
38 -30.82649004 26.91515923
39 26.31354313 -30.82649004
40 -13.69684737 26.31354313
41 0.90267055 -13.69684737
42 -3.82996813 0.90267055
43 15.88031058 -3.82996813
44 -1.34766344 15.88031058
45 -17.88031560 -1.34766344
46 -7.33315454 -17.88031560
47 13.65087545 -7.33315454
48 -27.98773831 13.65087545
49 22.32807601 -27.98773831
50 -7.03311489 22.32807601
51 14.70872358 -7.03311489
52 -29.31441343 14.70872358
53 -23.97108520 -29.31441343
54 -5.19901383 -23.97108520
55 1.97888868 -5.19901383
56 13.33954454 1.97888868
57 -6.18206542 13.33954454
58 -7.28885486 -6.18206542
59 -5.05796482 -7.28885486
60 NA -5.05796482
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -20.91692709 -10.79328689
[2,] -19.35415005 -20.91692709
[3,] -19.85882366 -19.35415005
[4,] -0.98449753 -19.85882366
[5,] 0.50024782 -0.98449753
[6,] -5.03015710 0.50024782
[7,] -11.63548077 -5.03015710
[8,] 18.39631464 -11.63548077
[9,] 16.61822724 18.39631464
[10,] -8.05433859 16.61822724
[11,] -12.03264222 -8.05433859
[12,] -10.14303087 -12.03264222
[13,] -16.58168833 -10.14303087
[14,] -9.99797019 -16.58168833
[15,] -16.41550723 -9.99797019
[16,] -11.29915586 -16.41550723
[17,] 6.94705283 -11.29915586
[18,] 4.74400625 6.94705283
[19,] 18.41592902 4.74400625
[20,] 26.20252558 18.41592902
[21,] 9.34742375 26.20252558
[22,] -0.04839387 9.34742375
[23,] 28.37078977 -0.04839387
[24,] 9.02986177 28.37078977
[25,] 34.74228468 9.02986177
[26,] -30.56865571 34.74228468
[27,] 32.98682251 -30.56865571
[28,] 25.10940370 32.98682251
[29,] 19.33127563 25.10940370
[30,] -4.15444901 19.33127563
[31,] -25.43532962 -4.15444901
[32,] 8.74399715 -25.43532962
[33,] 31.70883034 8.74399715
[34,] 0.66668186 31.70883034
[35,] -24.02090088 0.66668186
[36,] 16.39460905 -24.02090088
[37,] 26.91515923 16.39460905
[38,] -30.82649004 26.91515923
[39,] 26.31354313 -30.82649004
[40,] -13.69684737 26.31354313
[41,] 0.90267055 -13.69684737
[42,] -3.82996813 0.90267055
[43,] 15.88031058 -3.82996813
[44,] -1.34766344 15.88031058
[45,] -17.88031560 -1.34766344
[46,] -7.33315454 -17.88031560
[47,] 13.65087545 -7.33315454
[48,] -27.98773831 13.65087545
[49,] 22.32807601 -27.98773831
[50,] -7.03311489 22.32807601
[51,] 14.70872358 -7.03311489
[52,] -29.31441343 14.70872358
[53,] -23.97108520 -29.31441343
[54,] -5.19901383 -23.97108520
[55,] 1.97888868 -5.19901383
[56,] 13.33954454 1.97888868
[57,] -6.18206542 13.33954454
[58,] -7.28885486 -6.18206542
[59,] -5.05796482 -7.28885486
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -20.91692709 -10.79328689
2 -19.35415005 -20.91692709
3 -19.85882366 -19.35415005
4 -0.98449753 -19.85882366
5 0.50024782 -0.98449753
6 -5.03015710 0.50024782
7 -11.63548077 -5.03015710
8 18.39631464 -11.63548077
9 16.61822724 18.39631464
10 -8.05433859 16.61822724
11 -12.03264222 -8.05433859
12 -10.14303087 -12.03264222
13 -16.58168833 -10.14303087
14 -9.99797019 -16.58168833
15 -16.41550723 -9.99797019
16 -11.29915586 -16.41550723
17 6.94705283 -11.29915586
18 4.74400625 6.94705283
19 18.41592902 4.74400625
20 26.20252558 18.41592902
21 9.34742375 26.20252558
22 -0.04839387 9.34742375
23 28.37078977 -0.04839387
24 9.02986177 28.37078977
25 34.74228468 9.02986177
26 -30.56865571 34.74228468
27 32.98682251 -30.56865571
28 25.10940370 32.98682251
29 19.33127563 25.10940370
30 -4.15444901 19.33127563
31 -25.43532962 -4.15444901
32 8.74399715 -25.43532962
33 31.70883034 8.74399715
34 0.66668186 31.70883034
35 -24.02090088 0.66668186
36 16.39460905 -24.02090088
37 26.91515923 16.39460905
38 -30.82649004 26.91515923
39 26.31354313 -30.82649004
40 -13.69684737 26.31354313
41 0.90267055 -13.69684737
42 -3.82996813 0.90267055
43 15.88031058 -3.82996813
44 -1.34766344 15.88031058
45 -17.88031560 -1.34766344
46 -7.33315454 -17.88031560
47 13.65087545 -7.33315454
48 -27.98773831 13.65087545
49 22.32807601 -27.98773831
50 -7.03311489 22.32807601
51 14.70872358 -7.03311489
52 -29.31441343 14.70872358
53 -23.97108520 -29.31441343
54 -5.19901383 -23.97108520
55 1.97888868 -5.19901383
56 13.33954454 1.97888868
57 -6.18206542 13.33954454
58 -7.28885486 -6.18206542
59 -5.05796482 -7.28885486
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/7blzk1321784479.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/8t4qh1321784479.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/9t6od1321784479.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/rcomp/tmp/10cgta1321784479.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/11qi6b1321784479.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/12shhn1321784479.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/133otj1321784479.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/14esfp1321784479.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/rcomp/tmp/15efa71321784479.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/rcomp/tmp/16k1pd1321784479.tab")
+ }
>
> try(system("convert tmp/1nof71321784479.ps tmp/1nof71321784479.png",intern=TRUE))
character(0)
> try(system("convert tmp/2k6sf1321784479.ps tmp/2k6sf1321784479.png",intern=TRUE))
character(0)
> try(system("convert tmp/3x9pv1321784479.ps tmp/3x9pv1321784479.png",intern=TRUE))
character(0)
> try(system("convert tmp/4lab51321784479.ps tmp/4lab51321784479.png",intern=TRUE))
character(0)
> try(system("convert tmp/5vl651321784479.ps tmp/5vl651321784479.png",intern=TRUE))
character(0)
> try(system("convert tmp/6gpyv1321784479.ps tmp/6gpyv1321784479.png",intern=TRUE))
character(0)
> try(system("convert tmp/7blzk1321784479.ps tmp/7blzk1321784479.png",intern=TRUE))
character(0)
> try(system("convert tmp/8t4qh1321784479.ps tmp/8t4qh1321784479.png",intern=TRUE))
character(0)
> try(system("convert tmp/9t6od1321784479.ps tmp/9t6od1321784479.png",intern=TRUE))
character(0)
> try(system("convert tmp/10cgta1321784479.ps tmp/10cgta1321784479.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.040 0.300 4.298