R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
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(180144
+ ,966.2
+ ,173666
+ ,1153.2
+ ,165688
+ ,1328.3
+ ,161570
+ ,1144.5
+ ,156145
+ ,1477.1
+ ,153730
+ ,1234.9
+ ,182698
+ ,1119.1
+ ,200765
+ ,1356.9
+ ,176512
+ ,1217
+ ,166618
+ ,1440.5
+ ,158644
+ ,1556.6
+ ,159585
+ ,1303.6
+ ,163095
+ ,1421.5
+ ,159044
+ ,1172.5
+ ,155511
+ ,1422.1
+ ,153745
+ ,1263
+ ,150569
+ ,1428.1
+ ,150605
+ ,1347
+ ,179612
+ ,1224.2
+ ,194690
+ ,1201.3
+ ,189917
+ ,997.8
+ ,184128
+ ,1248.8
+ ,175335
+ ,1268.6
+ ,179566
+ ,1016.7
+ ,181140
+ ,1194.3
+ ,177876
+ ,1181.8
+ ,175041
+ ,1150.7
+ ,169292
+ ,1247.2
+ ,166070
+ ,1260.6
+ ,166972
+ ,1249.3
+ ,206348
+ ,1223.2
+ ,215706
+ ,1153
+ ,202108
+ ,1191.5
+ ,195411
+ ,1303.1
+ ,193111
+ ,1267.1
+ ,195198
+ ,1125.2
+ ,198770
+ ,1322.4
+ ,194163
+ ,1089.2
+ ,190420
+ ,1147.3
+ ,189733
+ ,1196.4
+ ,186029
+ ,1190.2
+ ,191531
+ ,1146
+ ,232571
+ ,1139.8
+ ,243477
+ ,1045.6
+ ,227247
+ ,1050.9
+ ,217859
+ ,1117.3
+ ,208679
+ ,1120
+ ,213188
+ ,1052.1
+ ,216234
+ ,1065.8
+ ,213586
+ ,1092.5
+ ,209465
+ ,1422
+ ,204045
+ ,1367.5
+ ,200237
+ ,1136.3
+ ,203666
+ ,1293.7
+ ,241476
+ ,1154.8
+ ,260307
+ ,1206.7
+ ,243324
+ ,1199
+ ,244460
+ ,1265
+ ,233575
+ ,1247.1
+ ,237217
+ ,1116.5
+ ,235243
+ ,1153.9
+ ,230354
+ ,1077.4
+ ,227184
+ ,1132.5
+ ,221678
+ ,1058.8
+ ,217142
+ ,1195.1
+ ,219452
+ ,1263.4
+ ,256446
+ ,1023.1
+ ,265845
+ ,1141
+ ,248624
+ ,1116.3
+ ,241114
+ ,1135.6
+ ,229245
+ ,1210.5
+ ,231805
+ ,1230
+ ,219277
+ ,1136.5
+ ,219313
+ ,1068.7
+ ,212610
+ ,1372.5
+ ,214771
+ ,1049.9
+ ,211142
+ ,1302.2
+ ,211457
+ ,1305.9
+ ,240048
+ ,1173.5
+ ,240636
+ ,1277.4
+ ,230580
+ ,1238.6
+ ,208795
+ ,1508.6
+ ,197922
+ ,1423.4
+ ,194596
+ ,1375.1)
+ ,dim=c(2
+ ,84)
+ ,dimnames=list(c('werkloosheid'
+ ,'Amerika')
+ ,1:84))
> y <- array(NA,dim=c(2,84),dimnames=list(c('werkloosheid','Amerika'),1:84))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> 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
werkloosheid Amerika M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 180144 966.2 1 0 0 0 0 0 0 0 0 0 0 1
2 173666 1153.2 0 1 0 0 0 0 0 0 0 0 0 2
3 165688 1328.3 0 0 1 0 0 0 0 0 0 0 0 3
4 161570 1144.5 0 0 0 1 0 0 0 0 0 0 0 4
5 156145 1477.1 0 0 0 0 1 0 0 0 0 0 0 5
6 153730 1234.9 0 0 0 0 0 1 0 0 0 0 0 6
7 182698 1119.1 0 0 0 0 0 0 1 0 0 0 0 7
8 200765 1356.9 0 0 0 0 0 0 0 1 0 0 0 8
9 176512 1217.0 0 0 0 0 0 0 0 0 1 0 0 9
10 166618 1440.5 0 0 0 0 0 0 0 0 0 1 0 10
11 158644 1556.6 0 0 0 0 0 0 0 0 0 0 1 11
12 159585 1303.6 0 0 0 0 0 0 0 0 0 0 0 12
13 163095 1421.5 1 0 0 0 0 0 0 0 0 0 0 13
14 159044 1172.5 0 1 0 0 0 0 0 0 0 0 0 14
15 155511 1422.1 0 0 1 0 0 0 0 0 0 0 0 15
16 153745 1263.0 0 0 0 1 0 0 0 0 0 0 0 16
17 150569 1428.1 0 0 0 0 1 0 0 0 0 0 0 17
18 150605 1347.0 0 0 0 0 0 1 0 0 0 0 0 18
19 179612 1224.2 0 0 0 0 0 0 1 0 0 0 0 19
20 194690 1201.3 0 0 0 0 0 0 0 1 0 0 0 20
21 189917 997.8 0 0 0 0 0 0 0 0 1 0 0 21
22 184128 1248.8 0 0 0 0 0 0 0 0 0 1 0 22
23 175335 1268.6 0 0 0 0 0 0 0 0 0 0 1 23
24 179566 1016.7 0 0 0 0 0 0 0 0 0 0 0 24
25 181140 1194.3 1 0 0 0 0 0 0 0 0 0 0 25
26 177876 1181.8 0 1 0 0 0 0 0 0 0 0 0 26
27 175041 1150.7 0 0 1 0 0 0 0 0 0 0 0 27
28 169292 1247.2 0 0 0 1 0 0 0 0 0 0 0 28
29 166070 1260.6 0 0 0 0 1 0 0 0 0 0 0 29
30 166972 1249.3 0 0 0 0 0 1 0 0 0 0 0 30
31 206348 1223.2 0 0 0 0 0 0 1 0 0 0 0 31
32 215706 1153.0 0 0 0 0 0 0 0 1 0 0 0 32
33 202108 1191.5 0 0 0 0 0 0 0 0 1 0 0 33
34 195411 1303.1 0 0 0 0 0 0 0 0 0 1 0 34
35 193111 1267.1 0 0 0 0 0 0 0 0 0 0 1 35
36 195198 1125.2 0 0 0 0 0 0 0 0 0 0 0 36
37 198770 1322.4 1 0 0 0 0 0 0 0 0 0 0 37
38 194163 1089.2 0 1 0 0 0 0 0 0 0 0 0 38
39 190420 1147.3 0 0 1 0 0 0 0 0 0 0 0 39
40 189733 1196.4 0 0 0 1 0 0 0 0 0 0 0 40
41 186029 1190.2 0 0 0 0 1 0 0 0 0 0 0 41
42 191531 1146.0 0 0 0 0 0 1 0 0 0 0 0 42
43 232571 1139.8 0 0 0 0 0 0 1 0 0 0 0 43
44 243477 1045.6 0 0 0 0 0 0 0 1 0 0 0 44
45 227247 1050.9 0 0 0 0 0 0 0 0 1 0 0 45
46 217859 1117.3 0 0 0 0 0 0 0 0 0 1 0 46
47 208679 1120.0 0 0 0 0 0 0 0 0 0 0 1 47
48 213188 1052.1 0 0 0 0 0 0 0 0 0 0 0 48
49 216234 1065.8 1 0 0 0 0 0 0 0 0 0 0 49
50 213586 1092.5 0 1 0 0 0 0 0 0 0 0 0 50
51 209465 1422.0 0 0 1 0 0 0 0 0 0 0 0 51
52 204045 1367.5 0 0 0 1 0 0 0 0 0 0 0 52
53 200237 1136.3 0 0 0 0 1 0 0 0 0 0 0 53
54 203666 1293.7 0 0 0 0 0 1 0 0 0 0 0 54
55 241476 1154.8 0 0 0 0 0 0 1 0 0 0 0 55
56 260307 1206.7 0 0 0 0 0 0 0 1 0 0 0 56
57 243324 1199.0 0 0 0 0 0 0 0 0 1 0 0 57
58 244460 1265.0 0 0 0 0 0 0 0 0 0 1 0 58
59 233575 1247.1 0 0 0 0 0 0 0 0 0 0 1 59
60 237217 1116.5 0 0 0 0 0 0 0 0 0 0 0 60
61 235243 1153.9 1 0 0 0 0 0 0 0 0 0 0 61
62 230354 1077.4 0 1 0 0 0 0 0 0 0 0 0 62
63 227184 1132.5 0 0 1 0 0 0 0 0 0 0 0 63
64 221678 1058.8 0 0 0 1 0 0 0 0 0 0 0 64
65 217142 1195.1 0 0 0 0 1 0 0 0 0 0 0 65
66 219452 1263.4 0 0 0 0 0 1 0 0 0 0 0 66
67 256446 1023.1 0 0 0 0 0 0 1 0 0 0 0 67
68 265845 1141.0 0 0 0 0 0 0 0 1 0 0 0 68
69 248624 1116.3 0 0 0 0 0 0 0 0 1 0 0 69
70 241114 1135.6 0 0 0 0 0 0 0 0 0 1 0 70
71 229245 1210.5 0 0 0 0 0 0 0 0 0 0 1 71
72 231805 1230.0 0 0 0 0 0 0 0 0 0 0 0 72
73 219277 1136.5 1 0 0 0 0 0 0 0 0 0 0 73
74 219313 1068.7 0 1 0 0 0 0 0 0 0 0 0 74
75 212610 1372.5 0 0 1 0 0 0 0 0 0 0 0 75
76 214771 1049.9 0 0 0 1 0 0 0 0 0 0 0 76
77 211142 1302.2 0 0 0 0 1 0 0 0 0 0 0 77
78 211457 1305.9 0 0 0 0 0 1 0 0 0 0 0 78
79 240048 1173.5 0 0 0 0 0 0 1 0 0 0 0 79
80 240636 1277.4 0 0 0 0 0 0 0 1 0 0 0 80
81 230580 1238.6 0 0 0 0 0 0 0 0 1 0 0 81
82 208795 1508.6 0 0 0 0 0 0 0 0 0 1 0 82
83 197922 1423.4 0 0 0 0 0 0 0 0 0 0 1 83
84 194596 1375.1 0 0 0 0 0 0 0 0 0 0 0 84
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Amerika M1 M2 M3 M4
234322.38 -63.99 7630.91 -840.21 4114.62 -5705.07
M5 M6 M7 M8 M9 M10
-4463.03 -5272.18 21232.59 35059.40 16056.04 15823.82
M11 t
6781.38 883.37
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-25942 -7615 -1238 7748 24021
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 234322.38 15365.74 15.250 < 2e-16 ***
Amerika -63.99 12.02 -5.324 1.17e-06 ***
M1 7630.91 6248.57 1.221 0.22610
M2 -840.21 6287.58 -0.134 0.89408
M3 4114.62 6358.74 0.647 0.51970
M4 -5705.07 6236.93 -0.915 0.36348
M5 -4463.03 6359.65 -0.702 0.48515
M6 -5272.18 6312.02 -0.835 0.40641
M7 21232.59 6235.11 3.405 0.00110 **
M8 35059.40 6229.87 5.628 3.51e-07 ***
M9 16056.04 6235.10 2.575 0.01214 *
M10 15823.82 6368.49 2.485 0.01536 *
M11 6781.38 6398.15 1.060 0.29284
t 883.37 53.53 16.501 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 11640 on 70 degrees of freedom
Multiple R-squared: 0.8683, Adjusted R-squared: 0.8438
F-statistic: 35.5 on 13 and 70 DF, p-value: < 2.2e-16
> 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.014279181 0.028558361 0.98572082
[2,] 0.010399376 0.020798753 0.98960062
[3,] 0.005237593 0.010475187 0.99476241
[4,] 0.001496168 0.002992336 0.99850383
[5,] 0.006858595 0.013717191 0.99314140
[6,] 0.012167617 0.024335234 0.98783238
[7,] 0.006597959 0.013195918 0.99340204
[8,] 0.004644797 0.009289594 0.99535520
[9,] 0.004839121 0.009678242 0.99516088
[10,] 0.005130131 0.010260263 0.99486987
[11,] 0.002987992 0.005975985 0.99701201
[12,] 0.003002235 0.006004469 0.99699777
[13,] 0.001490964 0.002981927 0.99850904
[14,] 0.001432706 0.002865412 0.99856729
[15,] 0.008577617 0.017155234 0.99142238
[16,] 0.006782976 0.013565951 0.99321702
[17,] 0.009662553 0.019325107 0.99033745
[18,] 0.008474836 0.016949671 0.99152516
[19,] 0.007351945 0.014703890 0.99264806
[20,] 0.009099579 0.018199158 0.99090042
[21,] 0.009473910 0.018947821 0.99052609
[22,] 0.007245413 0.014490827 0.99275459
[23,] 0.008490128 0.016980256 0.99150987
[24,] 0.008153679 0.016307357 0.99184632
[25,] 0.006935553 0.013871106 0.99306445
[26,] 0.011075470 0.022150940 0.98892453
[27,] 0.024900616 0.049801232 0.97509938
[28,] 0.031057229 0.062114459 0.96894277
[29,] 0.044236770 0.088473540 0.95576323
[30,] 0.067460969 0.134921937 0.93253903
[31,] 0.108112109 0.216224218 0.89188789
[32,] 0.223678145 0.447356290 0.77632186
[33,] 0.308782415 0.617564830 0.69121759
[34,] 0.338287260 0.676574521 0.66171274
[35,] 0.330484766 0.660969531 0.66951523
[36,] 0.306330124 0.612660248 0.69366988
[37,] 0.729213621 0.541572759 0.27078638
[38,] 0.889741203 0.220517593 0.11025880
[39,] 0.934042639 0.131914722 0.06595736
[40,] 0.932693377 0.134613246 0.06730662
[41,] 0.948480272 0.103039456 0.05151973
[42,] 0.939008178 0.121983644 0.06099182
[43,] 0.913351018 0.173297963 0.08664898
[44,] 0.887049813 0.225900374 0.11295019
[45,] 0.848058253 0.303883494 0.15194175
[46,] 0.772402988 0.455194023 0.22759701
[47,] 0.789044855 0.421910290 0.21095514
[48,] 0.685091853 0.629816293 0.31490815
[49,] 0.728972626 0.542054748 0.27102737
[50,] 0.720347263 0.559305475 0.27965274
[51,] 0.710404010 0.579191981 0.28959599
> postscript(file="/var/www/html/freestat/rcomp/tmp/1jq5m1229728413.ps",horizontal=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/html/freestat/rcomp/tmp/2pj0l1229728413.ps",horizontal=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/html/freestat/rcomp/tmp/32jrt1229728413.ps",horizontal=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/html/freestat/rcomp/tmp/460o41229728413.ps",horizontal=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/html/freestat/rcomp/tmp/589531229728413.ps",horizontal=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 = 84
Frequency = 1
1 2 3 4 5 6
-869.4330 12205.6856 9593.4280 2651.1383 16382.4513 -1604.1595
7 8 9 10 11 12
-7433.8745 11138.8062 -3945.8295 -190.1184 7423.7138 -1925.7143
13 14 15 16 17 18
613.9534 -11781.8014 -5182.1058 -8191.9426 -2929.2759 -8156.7505
19 20 21 22 23 24
-14395.3671 -15492.8231 -15166.9650 -5546.6404 -4913.6550 -10902.6986
25 26 27 28 29 30
-6479.0697 -2955.1478 -13618.3079 -4256.3363 -8746.3379 -8641.5934
31 32 33 34 35 36
1676.2312 -8167.7601 -1182.3026 -1389.6192 2165.9503 1071.3610
37 38 39 40 41 42
8747.1145 -3193.6624 -9057.2758 2333.7618 -3892.3644 -1292.7575
43 44 45 46 47 48
11962.3874 2130.7333 4359.8573 -1430.6241 -2278.7984 3783.5725
49 50 51 52 53 54
-808.0955 5840.0755 16964.2484 16993.3418 -3733.6228 9692.5513
55 56 57 58 59 60
11226.7608 18668.4538 19312.7605 24020.6847 20149.3999 21332.8518
61 62 63 64 65 66
13237.6508 11041.4720 5558.9007 4273.4638 6333.3352 12939.3613
67 68 69 70 71 72
7169.3957 9402.1612 8720.7068 1794.4873 2877.0984 12582.8412
73 74 75 76 77 78
-14442.1205 -11156.6215 -4258.8875 -13803.4268 -3414.1855 -2936.6517
79 80 81 82 83 84
-10205.5334 -17679.5712 -12098.2275 -17258.1698 -25423.7090 -25942.2135
> postscript(file="/var/www/html/freestat/rcomp/tmp/6zj661229728413.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 84
Frequency = 1
lag(myerror, k = 1) myerror
0 -869.4330 NA
1 12205.6856 -869.4330
2 9593.4280 12205.6856
3 2651.1383 9593.4280
4 16382.4513 2651.1383
5 -1604.1595 16382.4513
6 -7433.8745 -1604.1595
7 11138.8062 -7433.8745
8 -3945.8295 11138.8062
9 -190.1184 -3945.8295
10 7423.7138 -190.1184
11 -1925.7143 7423.7138
12 613.9534 -1925.7143
13 -11781.8014 613.9534
14 -5182.1058 -11781.8014
15 -8191.9426 -5182.1058
16 -2929.2759 -8191.9426
17 -8156.7505 -2929.2759
18 -14395.3671 -8156.7505
19 -15492.8231 -14395.3671
20 -15166.9650 -15492.8231
21 -5546.6404 -15166.9650
22 -4913.6550 -5546.6404
23 -10902.6986 -4913.6550
24 -6479.0697 -10902.6986
25 -2955.1478 -6479.0697
26 -13618.3079 -2955.1478
27 -4256.3363 -13618.3079
28 -8746.3379 -4256.3363
29 -8641.5934 -8746.3379
30 1676.2312 -8641.5934
31 -8167.7601 1676.2312
32 -1182.3026 -8167.7601
33 -1389.6192 -1182.3026
34 2165.9503 -1389.6192
35 1071.3610 2165.9503
36 8747.1145 1071.3610
37 -3193.6624 8747.1145
38 -9057.2758 -3193.6624
39 2333.7618 -9057.2758
40 -3892.3644 2333.7618
41 -1292.7575 -3892.3644
42 11962.3874 -1292.7575
43 2130.7333 11962.3874
44 4359.8573 2130.7333
45 -1430.6241 4359.8573
46 -2278.7984 -1430.6241
47 3783.5725 -2278.7984
48 -808.0955 3783.5725
49 5840.0755 -808.0955
50 16964.2484 5840.0755
51 16993.3418 16964.2484
52 -3733.6228 16993.3418
53 9692.5513 -3733.6228
54 11226.7608 9692.5513
55 18668.4538 11226.7608
56 19312.7605 18668.4538
57 24020.6847 19312.7605
58 20149.3999 24020.6847
59 21332.8518 20149.3999
60 13237.6508 21332.8518
61 11041.4720 13237.6508
62 5558.9007 11041.4720
63 4273.4638 5558.9007
64 6333.3352 4273.4638
65 12939.3613 6333.3352
66 7169.3957 12939.3613
67 9402.1612 7169.3957
68 8720.7068 9402.1612
69 1794.4873 8720.7068
70 2877.0984 1794.4873
71 12582.8412 2877.0984
72 -14442.1205 12582.8412
73 -11156.6215 -14442.1205
74 -4258.8875 -11156.6215
75 -13803.4268 -4258.8875
76 -3414.1855 -13803.4268
77 -2936.6517 -3414.1855
78 -10205.5334 -2936.6517
79 -17679.5712 -10205.5334
80 -12098.2275 -17679.5712
81 -17258.1698 -12098.2275
82 -25423.7090 -17258.1698
83 -25942.2135 -25423.7090
84 NA -25942.2135
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 12205.6856 -869.4330
[2,] 9593.4280 12205.6856
[3,] 2651.1383 9593.4280
[4,] 16382.4513 2651.1383
[5,] -1604.1595 16382.4513
[6,] -7433.8745 -1604.1595
[7,] 11138.8062 -7433.8745
[8,] -3945.8295 11138.8062
[9,] -190.1184 -3945.8295
[10,] 7423.7138 -190.1184
[11,] -1925.7143 7423.7138
[12,] 613.9534 -1925.7143
[13,] -11781.8014 613.9534
[14,] -5182.1058 -11781.8014
[15,] -8191.9426 -5182.1058
[16,] -2929.2759 -8191.9426
[17,] -8156.7505 -2929.2759
[18,] -14395.3671 -8156.7505
[19,] -15492.8231 -14395.3671
[20,] -15166.9650 -15492.8231
[21,] -5546.6404 -15166.9650
[22,] -4913.6550 -5546.6404
[23,] -10902.6986 -4913.6550
[24,] -6479.0697 -10902.6986
[25,] -2955.1478 -6479.0697
[26,] -13618.3079 -2955.1478
[27,] -4256.3363 -13618.3079
[28,] -8746.3379 -4256.3363
[29,] -8641.5934 -8746.3379
[30,] 1676.2312 -8641.5934
[31,] -8167.7601 1676.2312
[32,] -1182.3026 -8167.7601
[33,] -1389.6192 -1182.3026
[34,] 2165.9503 -1389.6192
[35,] 1071.3610 2165.9503
[36,] 8747.1145 1071.3610
[37,] -3193.6624 8747.1145
[38,] -9057.2758 -3193.6624
[39,] 2333.7618 -9057.2758
[40,] -3892.3644 2333.7618
[41,] -1292.7575 -3892.3644
[42,] 11962.3874 -1292.7575
[43,] 2130.7333 11962.3874
[44,] 4359.8573 2130.7333
[45,] -1430.6241 4359.8573
[46,] -2278.7984 -1430.6241
[47,] 3783.5725 -2278.7984
[48,] -808.0955 3783.5725
[49,] 5840.0755 -808.0955
[50,] 16964.2484 5840.0755
[51,] 16993.3418 16964.2484
[52,] -3733.6228 16993.3418
[53,] 9692.5513 -3733.6228
[54,] 11226.7608 9692.5513
[55,] 18668.4538 11226.7608
[56,] 19312.7605 18668.4538
[57,] 24020.6847 19312.7605
[58,] 20149.3999 24020.6847
[59,] 21332.8518 20149.3999
[60,] 13237.6508 21332.8518
[61,] 11041.4720 13237.6508
[62,] 5558.9007 11041.4720
[63,] 4273.4638 5558.9007
[64,] 6333.3352 4273.4638
[65,] 12939.3613 6333.3352
[66,] 7169.3957 12939.3613
[67,] 9402.1612 7169.3957
[68,] 8720.7068 9402.1612
[69,] 1794.4873 8720.7068
[70,] 2877.0984 1794.4873
[71,] 12582.8412 2877.0984
[72,] -14442.1205 12582.8412
[73,] -11156.6215 -14442.1205
[74,] -4258.8875 -11156.6215
[75,] -13803.4268 -4258.8875
[76,] -3414.1855 -13803.4268
[77,] -2936.6517 -3414.1855
[78,] -10205.5334 -2936.6517
[79,] -17679.5712 -10205.5334
[80,] -12098.2275 -17679.5712
[81,] -17258.1698 -12098.2275
[82,] -25423.7090 -17258.1698
[83,] -25942.2135 -25423.7090
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 12205.6856 -869.4330
2 9593.4280 12205.6856
3 2651.1383 9593.4280
4 16382.4513 2651.1383
5 -1604.1595 16382.4513
6 -7433.8745 -1604.1595
7 11138.8062 -7433.8745
8 -3945.8295 11138.8062
9 -190.1184 -3945.8295
10 7423.7138 -190.1184
11 -1925.7143 7423.7138
12 613.9534 -1925.7143
13 -11781.8014 613.9534
14 -5182.1058 -11781.8014
15 -8191.9426 -5182.1058
16 -2929.2759 -8191.9426
17 -8156.7505 -2929.2759
18 -14395.3671 -8156.7505
19 -15492.8231 -14395.3671
20 -15166.9650 -15492.8231
21 -5546.6404 -15166.9650
22 -4913.6550 -5546.6404
23 -10902.6986 -4913.6550
24 -6479.0697 -10902.6986
25 -2955.1478 -6479.0697
26 -13618.3079 -2955.1478
27 -4256.3363 -13618.3079
28 -8746.3379 -4256.3363
29 -8641.5934 -8746.3379
30 1676.2312 -8641.5934
31 -8167.7601 1676.2312
32 -1182.3026 -8167.7601
33 -1389.6192 -1182.3026
34 2165.9503 -1389.6192
35 1071.3610 2165.9503
36 8747.1145 1071.3610
37 -3193.6624 8747.1145
38 -9057.2758 -3193.6624
39 2333.7618 -9057.2758
40 -3892.3644 2333.7618
41 -1292.7575 -3892.3644
42 11962.3874 -1292.7575
43 2130.7333 11962.3874
44 4359.8573 2130.7333
45 -1430.6241 4359.8573
46 -2278.7984 -1430.6241
47 3783.5725 -2278.7984
48 -808.0955 3783.5725
49 5840.0755 -808.0955
50 16964.2484 5840.0755
51 16993.3418 16964.2484
52 -3733.6228 16993.3418
53 9692.5513 -3733.6228
54 11226.7608 9692.5513
55 18668.4538 11226.7608
56 19312.7605 18668.4538
57 24020.6847 19312.7605
58 20149.3999 24020.6847
59 21332.8518 20149.3999
60 13237.6508 21332.8518
61 11041.4720 13237.6508
62 5558.9007 11041.4720
63 4273.4638 5558.9007
64 6333.3352 4273.4638
65 12939.3613 6333.3352
66 7169.3957 12939.3613
67 9402.1612 7169.3957
68 8720.7068 9402.1612
69 1794.4873 8720.7068
70 2877.0984 1794.4873
71 12582.8412 2877.0984
72 -14442.1205 12582.8412
73 -11156.6215 -14442.1205
74 -4258.8875 -11156.6215
75 -13803.4268 -4258.8875
76 -3414.1855 -13803.4268
77 -2936.6517 -3414.1855
78 -10205.5334 -2936.6517
79 -17679.5712 -10205.5334
80 -12098.2275 -17679.5712
81 -17258.1698 -12098.2275
82 -25423.7090 -17258.1698
83 -25942.2135 -25423.7090
> 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/html/freestat/rcomp/tmp/7v2up1229728413.ps",horizontal=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/html/freestat/rcomp/tmp/8lzm51229728413.ps",horizontal=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/html/freestat/rcomp/tmp/9d3181229728413.ps",horizontal=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/html/freestat/rcomp/tmp/1017g61229728413.ps",horizontal=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/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/html/freestat/rcomp/tmp/11t7or1229728414.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/html/freestat/rcomp/tmp/12wd7j1229728414.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/html/freestat/rcomp/tmp/13mf4y1229728414.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/html/freestat/rcomp/tmp/14m1tr1229728414.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/html/freestat/rcomp/tmp/15ms1m1229728414.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/html/freestat/rcomp/tmp/162r621229728414.tab")
+ }
>
> system("convert tmp/1jq5m1229728413.ps tmp/1jq5m1229728413.png")
> system("convert tmp/2pj0l1229728413.ps tmp/2pj0l1229728413.png")
> system("convert tmp/32jrt1229728413.ps tmp/32jrt1229728413.png")
> system("convert tmp/460o41229728413.ps tmp/460o41229728413.png")
> system("convert tmp/589531229728413.ps tmp/589531229728413.png")
> system("convert tmp/6zj661229728413.ps tmp/6zj661229728413.png")
> system("convert tmp/7v2up1229728413.ps tmp/7v2up1229728413.png")
> system("convert tmp/8lzm51229728413.ps tmp/8lzm51229728413.png")
> system("convert tmp/9d3181229728413.ps tmp/9d3181229728413.png")
> system("convert tmp/1017g61229728413.ps tmp/1017g61229728413.png")
>
>
> proc.time()
user system elapsed
4.099 2.556 7.674