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 = 'No Linear Trend'
> par2 = 'Do not include Seasonal 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
1 180144 966.2
2 173666 1153.2
3 165688 1328.3
4 161570 1144.5
5 156145 1477.1
6 153730 1234.9
7 182698 1119.1
8 200765 1356.9
9 176512 1217.0
10 166618 1440.5
11 158644 1556.6
12 159585 1303.6
13 163095 1421.5
14 159044 1172.5
15 155511 1422.1
16 153745 1263.0
17 150569 1428.1
18 150605 1347.0
19 179612 1224.2
20 194690 1201.3
21 189917 997.8
22 184128 1248.8
23 175335 1268.6
24 179566 1016.7
25 181140 1194.3
26 177876 1181.8
27 175041 1150.7
28 169292 1247.2
29 166070 1260.6
30 166972 1249.3
31 206348 1223.2
32 215706 1153.0
33 202108 1191.5
34 195411 1303.1
35 193111 1267.1
36 195198 1125.2
37 198770 1322.4
38 194163 1089.2
39 190420 1147.3
40 189733 1196.4
41 186029 1190.2
42 191531 1146.0
43 232571 1139.8
44 243477 1045.6
45 227247 1050.9
46 217859 1117.3
47 208679 1120.0
48 213188 1052.1
49 216234 1065.8
50 213586 1092.5
51 209465 1422.0
52 204045 1367.5
53 200237 1136.3
54 203666 1293.7
55 241476 1154.8
56 260307 1206.7
57 243324 1199.0
58 244460 1265.0
59 233575 1247.1
60 237217 1116.5
61 235243 1153.9
62 230354 1077.4
63 227184 1132.5
64 221678 1058.8
65 217142 1195.1
66 219452 1263.4
67 256446 1023.1
68 265845 1141.0
69 248624 1116.3
70 241114 1135.6
71 229245 1210.5
72 231805 1230.0
73 219277 1136.5
74 219313 1068.7
75 212610 1372.5
76 214771 1049.9
77 211142 1302.2
78 211457 1305.9
79 240048 1173.5
80 240636 1277.4
81 230580 1238.6
82 208795 1508.6
83 197922 1423.4
84 194596 1375.1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Amerika
314745.6 -93.1
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-46633.2 -21593.5 623.3 22720.5 57894.0
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 314745.61 29367.57 10.717 < 2e-16 ***
Amerika -93.09 24.06 -3.869 0.000218 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 27250 on 82 degrees of freedom
Multiple R-squared: 0.1544, Adjusted R-squared: 0.1441
F-statistic: 14.97 on 1 and 82 DF, p-value: 0.000218
> 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.014890663 0.029781326 0.985109337
[2,] 0.014388380 0.028776760 0.985611620
[3,] 0.011163473 0.022326945 0.988836527
[4,] 0.128600674 0.257201348 0.871399326
[5,] 0.074660852 0.149321704 0.925339148
[6,] 0.039227082 0.078454164 0.960772918
[7,] 0.020573186 0.041146373 0.979426814
[8,] 0.013282861 0.026565722 0.986717139
[9,] 0.006610111 0.013220222 0.993389889
[10,] 0.005739137 0.011478274 0.994260863
[11,] 0.003664014 0.007328029 0.996335986
[12,] 0.003763705 0.007527409 0.996236295
[13,] 0.003223657 0.006447314 0.996776343
[14,] 0.003522748 0.007045495 0.996477252
[15,] 0.002831344 0.005662688 0.997168656
[16,] 0.005334429 0.010668858 0.994665571
[17,] 0.004131053 0.008262107 0.995868947
[18,] 0.003680875 0.007361749 0.996319125
[19,] 0.002698306 0.005396612 0.997301694
[20,] 0.002293026 0.004586053 0.997706974
[21,] 0.001853713 0.003707426 0.998146287
[22,] 0.001523083 0.003046165 0.998476917
[23,] 0.001450863 0.002901726 0.998549137
[24,] 0.001497158 0.002994317 0.998502842
[25,] 0.001942268 0.003884536 0.998057732
[26,] 0.002913488 0.005826976 0.997086512
[27,] 0.013893370 0.027786739 0.986106630
[28,] 0.050831472 0.101662944 0.949168528
[29,] 0.069772142 0.139544283 0.930227858
[30,] 0.089082346 0.178164693 0.910917654
[31,] 0.101301660 0.202603320 0.898698340
[32,] 0.109848047 0.219696094 0.890151953
[33,] 0.140701942 0.281403884 0.859298058
[34,] 0.158300040 0.316600080 0.841699960
[35,] 0.187579649 0.375159298 0.812420351
[36,] 0.225033975 0.450067950 0.774966025
[37,] 0.298950950 0.597901901 0.701049050
[38,] 0.388732905 0.777465809 0.611267095
[39,] 0.600868377 0.798263247 0.399131623
[40,] 0.766370359 0.467259281 0.233629641
[41,] 0.782913721 0.434172558 0.217086279
[42,] 0.788805282 0.422389437 0.211194718
[43,] 0.800337898 0.399324204 0.199662102
[44,] 0.812620893 0.374758214 0.187379107
[45,] 0.821894717 0.356210566 0.178105283
[46,] 0.838823805 0.322352391 0.161176195
[47,] 0.877675806 0.244648388 0.122324194
[48,] 0.884444323 0.231111354 0.115555677
[49,] 0.930745993 0.138508014 0.069254007
[50,] 0.937974225 0.124051550 0.062025775
[51,] 0.955560172 0.088879657 0.044439828
[52,] 0.994021288 0.011957424 0.005978712
[53,] 0.996257051 0.007485898 0.003742949
[54,] 0.998572848 0.002854304 0.001427152
[55,] 0.998562001 0.002875997 0.001437999
[56,] 0.997972804 0.004054393 0.002027196
[57,] 0.997176725 0.005646551 0.002823275
[58,] 0.995545689 0.008908621 0.004454311
[59,] 0.993018989 0.013962022 0.006981011
[60,] 0.992275249 0.015449502 0.007724751
[61,] 0.989277803 0.021444395 0.010722197
[62,] 0.983434393 0.033131215 0.016565607
[63,] 0.980523566 0.038952868 0.019476434
[64,] 0.996518970 0.006962060 0.003481030
[65,] 0.997063329 0.005873341 0.002936671
[66,] 0.996477807 0.007044386 0.003522193
[67,] 0.993828845 0.012342311 0.006171155
[68,] 0.991329718 0.017340564 0.008670282
[69,] 0.982125634 0.035748732 0.017874366
[70,] 0.969095881 0.061808238 0.030904119
[71,] 0.941296964 0.117406072 0.058703036
[72,] 0.965528541 0.068942918 0.034471459
[73,] 0.936930158 0.126139684 0.063069842
[74,] 0.888265591 0.223468817 0.111734409
[75,] 0.772568059 0.454863883 0.227431941
> postscript(file="/var/www/html/rcomp/tmp/1pe6x1229727293.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/rcomp/tmp/2m3hs1229727293.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/rcomp/tmp/3arui1229727293.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/rcomp/tmp/4em841229727293.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/rcomp/tmp/5qjct1229727293.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
-44657.2864 -33727.3086 -25405.1111 -46633.1984 -21096.2014 -46057.7910
7 8 9 10 11 12
-27869.7045 12334.2855 -24942.1161 -14030.3244 -11196.4836 -33807.4536
13 14 15 16 17 18
-19322.0494 -46552.6563 -26850.1949 -43426.9397 -31233.6502 -38747.3133
19 20 21 22 23 24
-21171.8624 -8225.6415 -31942.6174 -14365.8290 -21315.6313 -40534.2015
25 26 27 28 29 30
-22427.2771 -26854.9119 -32585.0355 -29350.7742 -31325.3576 -31475.2836
31 32 33 34 35 36
5471.0468 8294.0733 -1719.9313 1972.0010 -3679.2675 -14801.8507
37 38 39 40 41 42
7127.6532 -19188.1192 -17522.5442 -13638.7864 -17919.9493 -16532.5622
43 44 45 46 47 48
23930.2748 26067.1224 10330.5036 7123.7321 -1804.9228 -3616.7875
49 50 51 52 53 54
704.5563 542.0804 27094.4960 16601.0479 -8729.5429 9351.9475
55 56 57 58 59 60
34231.6367 57894.0487 40194.2497 47474.2418 34922.9167 26407.2594
61 62 63 64 65 66
27914.8550 15904.4095 17863.7121 5496.9208 13649.1956 22317.2966
67 68 69 70 71 72
36941.5796 57315.9838 37795.6413 32082.2935 27185.7938 31561.0642
73 74 75 76 77 78
10329.0752 4053.5196 25631.5019 -2238.5872 17619.2193 18278.6552
79 80 81 82 83 84
34544.4345 44804.5677 31136.6450 34486.1585 15681.8231 7859.5379
> postscript(file="/var/www/html/rcomp/tmp/6u7w61229727293.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 -44657.2864 NA
1 -33727.3086 -44657.2864
2 -25405.1111 -33727.3086
3 -46633.1984 -25405.1111
4 -21096.2014 -46633.1984
5 -46057.7910 -21096.2014
6 -27869.7045 -46057.7910
7 12334.2855 -27869.7045
8 -24942.1161 12334.2855
9 -14030.3244 -24942.1161
10 -11196.4836 -14030.3244
11 -33807.4536 -11196.4836
12 -19322.0494 -33807.4536
13 -46552.6563 -19322.0494
14 -26850.1949 -46552.6563
15 -43426.9397 -26850.1949
16 -31233.6502 -43426.9397
17 -38747.3133 -31233.6502
18 -21171.8624 -38747.3133
19 -8225.6415 -21171.8624
20 -31942.6174 -8225.6415
21 -14365.8290 -31942.6174
22 -21315.6313 -14365.8290
23 -40534.2015 -21315.6313
24 -22427.2771 -40534.2015
25 -26854.9119 -22427.2771
26 -32585.0355 -26854.9119
27 -29350.7742 -32585.0355
28 -31325.3576 -29350.7742
29 -31475.2836 -31325.3576
30 5471.0468 -31475.2836
31 8294.0733 5471.0468
32 -1719.9313 8294.0733
33 1972.0010 -1719.9313
34 -3679.2675 1972.0010
35 -14801.8507 -3679.2675
36 7127.6532 -14801.8507
37 -19188.1192 7127.6532
38 -17522.5442 -19188.1192
39 -13638.7864 -17522.5442
40 -17919.9493 -13638.7864
41 -16532.5622 -17919.9493
42 23930.2748 -16532.5622
43 26067.1224 23930.2748
44 10330.5036 26067.1224
45 7123.7321 10330.5036
46 -1804.9228 7123.7321
47 -3616.7875 -1804.9228
48 704.5563 -3616.7875
49 542.0804 704.5563
50 27094.4960 542.0804
51 16601.0479 27094.4960
52 -8729.5429 16601.0479
53 9351.9475 -8729.5429
54 34231.6367 9351.9475
55 57894.0487 34231.6367
56 40194.2497 57894.0487
57 47474.2418 40194.2497
58 34922.9167 47474.2418
59 26407.2594 34922.9167
60 27914.8550 26407.2594
61 15904.4095 27914.8550
62 17863.7121 15904.4095
63 5496.9208 17863.7121
64 13649.1956 5496.9208
65 22317.2966 13649.1956
66 36941.5796 22317.2966
67 57315.9838 36941.5796
68 37795.6413 57315.9838
69 32082.2935 37795.6413
70 27185.7938 32082.2935
71 31561.0642 27185.7938
72 10329.0752 31561.0642
73 4053.5196 10329.0752
74 25631.5019 4053.5196
75 -2238.5872 25631.5019
76 17619.2193 -2238.5872
77 18278.6552 17619.2193
78 34544.4345 18278.6552
79 44804.5677 34544.4345
80 31136.6450 44804.5677
81 34486.1585 31136.6450
82 15681.8231 34486.1585
83 7859.5379 15681.8231
84 NA 7859.5379
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -33727.3086 -44657.2864
[2,] -25405.1111 -33727.3086
[3,] -46633.1984 -25405.1111
[4,] -21096.2014 -46633.1984
[5,] -46057.7910 -21096.2014
[6,] -27869.7045 -46057.7910
[7,] 12334.2855 -27869.7045
[8,] -24942.1161 12334.2855
[9,] -14030.3244 -24942.1161
[10,] -11196.4836 -14030.3244
[11,] -33807.4536 -11196.4836
[12,] -19322.0494 -33807.4536
[13,] -46552.6563 -19322.0494
[14,] -26850.1949 -46552.6563
[15,] -43426.9397 -26850.1949
[16,] -31233.6502 -43426.9397
[17,] -38747.3133 -31233.6502
[18,] -21171.8624 -38747.3133
[19,] -8225.6415 -21171.8624
[20,] -31942.6174 -8225.6415
[21,] -14365.8290 -31942.6174
[22,] -21315.6313 -14365.8290
[23,] -40534.2015 -21315.6313
[24,] -22427.2771 -40534.2015
[25,] -26854.9119 -22427.2771
[26,] -32585.0355 -26854.9119
[27,] -29350.7742 -32585.0355
[28,] -31325.3576 -29350.7742
[29,] -31475.2836 -31325.3576
[30,] 5471.0468 -31475.2836
[31,] 8294.0733 5471.0468
[32,] -1719.9313 8294.0733
[33,] 1972.0010 -1719.9313
[34,] -3679.2675 1972.0010
[35,] -14801.8507 -3679.2675
[36,] 7127.6532 -14801.8507
[37,] -19188.1192 7127.6532
[38,] -17522.5442 -19188.1192
[39,] -13638.7864 -17522.5442
[40,] -17919.9493 -13638.7864
[41,] -16532.5622 -17919.9493
[42,] 23930.2748 -16532.5622
[43,] 26067.1224 23930.2748
[44,] 10330.5036 26067.1224
[45,] 7123.7321 10330.5036
[46,] -1804.9228 7123.7321
[47,] -3616.7875 -1804.9228
[48,] 704.5563 -3616.7875
[49,] 542.0804 704.5563
[50,] 27094.4960 542.0804
[51,] 16601.0479 27094.4960
[52,] -8729.5429 16601.0479
[53,] 9351.9475 -8729.5429
[54,] 34231.6367 9351.9475
[55,] 57894.0487 34231.6367
[56,] 40194.2497 57894.0487
[57,] 47474.2418 40194.2497
[58,] 34922.9167 47474.2418
[59,] 26407.2594 34922.9167
[60,] 27914.8550 26407.2594
[61,] 15904.4095 27914.8550
[62,] 17863.7121 15904.4095
[63,] 5496.9208 17863.7121
[64,] 13649.1956 5496.9208
[65,] 22317.2966 13649.1956
[66,] 36941.5796 22317.2966
[67,] 57315.9838 36941.5796
[68,] 37795.6413 57315.9838
[69,] 32082.2935 37795.6413
[70,] 27185.7938 32082.2935
[71,] 31561.0642 27185.7938
[72,] 10329.0752 31561.0642
[73,] 4053.5196 10329.0752
[74,] 25631.5019 4053.5196
[75,] -2238.5872 25631.5019
[76,] 17619.2193 -2238.5872
[77,] 18278.6552 17619.2193
[78,] 34544.4345 18278.6552
[79,] 44804.5677 34544.4345
[80,] 31136.6450 44804.5677
[81,] 34486.1585 31136.6450
[82,] 15681.8231 34486.1585
[83,] 7859.5379 15681.8231
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -33727.3086 -44657.2864
2 -25405.1111 -33727.3086
3 -46633.1984 -25405.1111
4 -21096.2014 -46633.1984
5 -46057.7910 -21096.2014
6 -27869.7045 -46057.7910
7 12334.2855 -27869.7045
8 -24942.1161 12334.2855
9 -14030.3244 -24942.1161
10 -11196.4836 -14030.3244
11 -33807.4536 -11196.4836
12 -19322.0494 -33807.4536
13 -46552.6563 -19322.0494
14 -26850.1949 -46552.6563
15 -43426.9397 -26850.1949
16 -31233.6502 -43426.9397
17 -38747.3133 -31233.6502
18 -21171.8624 -38747.3133
19 -8225.6415 -21171.8624
20 -31942.6174 -8225.6415
21 -14365.8290 -31942.6174
22 -21315.6313 -14365.8290
23 -40534.2015 -21315.6313
24 -22427.2771 -40534.2015
25 -26854.9119 -22427.2771
26 -32585.0355 -26854.9119
27 -29350.7742 -32585.0355
28 -31325.3576 -29350.7742
29 -31475.2836 -31325.3576
30 5471.0468 -31475.2836
31 8294.0733 5471.0468
32 -1719.9313 8294.0733
33 1972.0010 -1719.9313
34 -3679.2675 1972.0010
35 -14801.8507 -3679.2675
36 7127.6532 -14801.8507
37 -19188.1192 7127.6532
38 -17522.5442 -19188.1192
39 -13638.7864 -17522.5442
40 -17919.9493 -13638.7864
41 -16532.5622 -17919.9493
42 23930.2748 -16532.5622
43 26067.1224 23930.2748
44 10330.5036 26067.1224
45 7123.7321 10330.5036
46 -1804.9228 7123.7321
47 -3616.7875 -1804.9228
48 704.5563 -3616.7875
49 542.0804 704.5563
50 27094.4960 542.0804
51 16601.0479 27094.4960
52 -8729.5429 16601.0479
53 9351.9475 -8729.5429
54 34231.6367 9351.9475
55 57894.0487 34231.6367
56 40194.2497 57894.0487
57 47474.2418 40194.2497
58 34922.9167 47474.2418
59 26407.2594 34922.9167
60 27914.8550 26407.2594
61 15904.4095 27914.8550
62 17863.7121 15904.4095
63 5496.9208 17863.7121
64 13649.1956 5496.9208
65 22317.2966 13649.1956
66 36941.5796 22317.2966
67 57315.9838 36941.5796
68 37795.6413 57315.9838
69 32082.2935 37795.6413
70 27185.7938 32082.2935
71 31561.0642 27185.7938
72 10329.0752 31561.0642
73 4053.5196 10329.0752
74 25631.5019 4053.5196
75 -2238.5872 25631.5019
76 17619.2193 -2238.5872
77 18278.6552 17619.2193
78 34544.4345 18278.6552
79 44804.5677 34544.4345
80 31136.6450 44804.5677
81 34486.1585 31136.6450
82 15681.8231 34486.1585
83 7859.5379 15681.8231
> 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/rcomp/tmp/7iaog1229727293.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/rcomp/tmp/8bm5v1229727293.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/rcomp/tmp/98yh81229727293.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/rcomp/tmp/10zzjk1229727293.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/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/rcomp/tmp/11919f1229727293.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/rcomp/tmp/12dci81229727293.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/rcomp/tmp/133nyv1229727293.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/rcomp/tmp/143hdy1229727293.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/rcomp/tmp/15xwgc1229727293.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/rcomp/tmp/16qe3r1229727293.tab")
+ }
>
> system("convert tmp/1pe6x1229727293.ps tmp/1pe6x1229727293.png")
> system("convert tmp/2m3hs1229727293.ps tmp/2m3hs1229727293.png")
> system("convert tmp/3arui1229727293.ps tmp/3arui1229727293.png")
> system("convert tmp/4em841229727293.ps tmp/4em841229727293.png")
> system("convert tmp/5qjct1229727293.ps tmp/5qjct1229727293.png")
> system("convert tmp/6u7w61229727293.ps tmp/6u7w61229727293.png")
> system("convert tmp/7iaog1229727293.ps tmp/7iaog1229727293.png")
> system("convert tmp/8bm5v1229727293.ps tmp/8bm5v1229727293.png")
> system("convert tmp/98yh81229727293.ps tmp/98yh81229727293.png")
> system("convert tmp/10zzjk1229727293.ps tmp/10zzjk1229727293.png")
>
>
> proc.time()
user system elapsed
2.821 1.669 4.782