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
+ ,11554.5
+ ,173666
+ ,13182.1
+ ,165688
+ ,14800.1
+ ,161570
+ ,12150.7
+ ,156145
+ ,14478.2
+ ,153730
+ ,13253.9
+ ,182698
+ ,12036.8
+ ,200765
+ ,12653.2
+ ,176512
+ ,14035.4
+ ,166618
+ ,14571.4
+ ,158644
+ ,15400.9
+ ,159585
+ ,14283.2
+ ,163095
+ ,14485.3
+ ,159044
+ ,14196.3
+ ,155511
+ ,15559.1
+ ,153745
+ ,13767.4
+ ,150569
+ ,14634
+ ,150605
+ ,14381.1
+ ,179612
+ ,12509.9
+ ,194690
+ ,12122.3
+ ,189917
+ ,13122.3
+ ,184128
+ ,13908.7
+ ,175335
+ ,13456.5
+ ,179566
+ ,12441.6
+ ,181140
+ ,12953
+ ,177876
+ ,13057.2
+ ,175041
+ ,14350.1
+ ,169292
+ ,13830.2
+ ,166070
+ ,13755.5
+ ,166972
+ ,13574.4
+ ,206348
+ ,12802.6
+ ,215706
+ ,11737.3
+ ,202108
+ ,13850.2
+ ,195411
+ ,15081.8
+ ,193111
+ ,13653.3
+ ,195198
+ ,14019.1
+ ,198770
+ ,13962
+ ,194163
+ ,13768.7
+ ,190420
+ ,14747.1
+ ,189733
+ ,13858.1
+ ,186029
+ ,13188
+ ,191531
+ ,13693.1
+ ,232571
+ ,12970
+ ,243477
+ ,11392.8
+ ,227247
+ ,13985.2
+ ,217859
+ ,14994.7
+ ,208679
+ ,13584.7
+ ,213188
+ ,14257.8
+ ,216234
+ ,13553.4
+ ,213586
+ ,14007.3
+ ,209465
+ ,16535.8
+ ,204045
+ ,14721.4
+ ,200237
+ ,13664.6
+ ,203666
+ ,16405.9
+ ,241476
+ ,13829.4
+ ,260307
+ ,13735.6
+ ,243324
+ ,15870.5
+ ,244460
+ ,15962.4
+ ,233575
+ ,15744.1
+ ,237217
+ ,16083.7
+ ,235243
+ ,14863.9
+ ,230354
+ ,15533.1
+ ,227184
+ ,17473.1
+ ,221678
+ ,15925.5
+ ,217142
+ ,15573.7
+ ,219452
+ ,17495
+ ,256446
+ ,14155.8
+ ,265845
+ ,14913.9
+ ,248624
+ ,17250.4
+ ,241114
+ ,15879.8
+ ,229245
+ ,17647.8
+ ,231805
+ ,17749.9
+ ,219277
+ ,17111.8
+ ,219313
+ ,16934.8
+ ,212610
+ ,20280
+ ,214771
+ ,16238.2
+ ,211142
+ ,17896.1
+ ,211457
+ ,18089.3
+ ,240048
+ ,15660
+ ,240636
+ ,16162.4
+ ,230580
+ ,17850.1
+ ,208795
+ ,18520.4
+ ,197922
+ ,18524.7
+ ,194596
+ ,16843.7)
+ ,dim=c(2
+ ,84)
+ ,dimnames=list(c('werkloosheid'
+ ,'invoer')
+ ,1:84))
> y <- array(NA,dim=c(2,84),dimnames=list(c('werkloosheid','invoer'),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 invoer
1 180144 11554.5
2 173666 13182.1
3 165688 14800.1
4 161570 12150.7
5 156145 14478.2
6 153730 13253.9
7 182698 12036.8
8 200765 12653.2
9 176512 14035.4
10 166618 14571.4
11 158644 15400.9
12 159585 14283.2
13 163095 14485.3
14 159044 14196.3
15 155511 15559.1
16 153745 13767.4
17 150569 14634.0
18 150605 14381.1
19 179612 12509.9
20 194690 12122.3
21 189917 13122.3
22 184128 13908.7
23 175335 13456.5
24 179566 12441.6
25 181140 12953.0
26 177876 13057.2
27 175041 14350.1
28 169292 13830.2
29 166070 13755.5
30 166972 13574.4
31 206348 12802.6
32 215706 11737.3
33 202108 13850.2
34 195411 15081.8
35 193111 13653.3
36 195198 14019.1
37 198770 13962.0
38 194163 13768.7
39 190420 14747.1
40 189733 13858.1
41 186029 13188.0
42 191531 13693.1
43 232571 12970.0
44 243477 11392.8
45 227247 13985.2
46 217859 14994.7
47 208679 13584.7
48 213188 14257.8
49 216234 13553.4
50 213586 14007.3
51 209465 16535.8
52 204045 14721.4
53 200237 13664.6
54 203666 16405.9
55 241476 13829.4
56 260307 13735.6
57 243324 15870.5
58 244460 15962.4
59 233575 15744.1
60 237217 16083.7
61 235243 14863.9
62 230354 15533.1
63 227184 17473.1
64 221678 15925.5
65 217142 15573.7
66 219452 17495.0
67 256446 14155.8
68 265845 14913.9
69 248624 17250.4
70 241114 15879.8
71 229245 17647.8
72 231805 17749.9
73 219277 17111.8
74 219313 16934.8
75 212610 20280.0
76 214771 16238.2
77 211142 17896.1
78 211457 18089.3
79 240048 15660.0
80 240636 16162.4
81 230580 17850.1
82 208795 18520.4
83 197922 18524.7
84 194596 16843.7
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) invoer
1.175e+05 5.709e+00
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-50821 -19076 -2342 15000 64385
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.175e+05 2.481e+04 4.736 9e-06 ***
invoer 5.709e+00 1.670e+00 3.419 0.000981 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 27720 on 82 degrees of freedom
Multiple R-squared: 0.1248, Adjusted R-squared: 0.1141
F-statistic: 11.69 on 1 and 82 DF, p-value: 0.0009815
> 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.0423894702 0.0847789405 0.9576105298
[2,] 0.0294500809 0.0589001618 0.9705499191
[3,] 0.0163299004 0.0326598008 0.9836700996
[4,] 0.0601427737 0.1202855473 0.9398572263
[5,] 0.0355329012 0.0710658024 0.9644670988
[6,] 0.0176070689 0.0352141379 0.9823929311
[7,] 0.0090751559 0.0181503117 0.9909248441
[8,] 0.0049283587 0.0098567173 0.9950716413
[9,] 0.0023951806 0.0047903611 0.9976048194
[10,] 0.0013587537 0.0027175074 0.9986412463
[11,] 0.0007879264 0.0015758528 0.9992120736
[12,] 0.0008325253 0.0016650506 0.9991674747
[13,] 0.0008324842 0.0016649684 0.9991675158
[14,] 0.0009685320 0.0019370641 0.9990314680
[15,] 0.0005178536 0.0010357072 0.9994821464
[16,] 0.0005237262 0.0010474524 0.9994762738
[17,] 0.0006229563 0.0012459126 0.9993770437
[18,] 0.0007624406 0.0015248811 0.9992375594
[19,] 0.0004926556 0.0009853112 0.9995073444
[20,] 0.0002772874 0.0005545748 0.9997227126
[21,] 0.0001737186 0.0003474372 0.9998262814
[22,] 0.0001113443 0.0002226885 0.9998886557
[23,] 0.0001193076 0.0002386152 0.9998806924
[24,] 0.0001112739 0.0002225479 0.9998887261
[25,] 0.0001379355 0.0002758710 0.9998620645
[26,] 0.0002025833 0.0004051665 0.9997974167
[27,] 0.0009592051 0.0019184101 0.9990407949
[28,] 0.0020967615 0.0041935231 0.9979032385
[29,] 0.0073646737 0.0147293475 0.9926353263
[30,] 0.0263898004 0.0527796008 0.9736101996
[31,] 0.0325488544 0.0650977088 0.9674511456
[32,] 0.0461907355 0.0923814711 0.9538092645
[33,] 0.0658118469 0.1316236937 0.9341881531
[34,] 0.0815269305 0.1630538609 0.9184730695
[35,] 0.1186810507 0.2373621014 0.8813189493
[36,] 0.1557479537 0.3114959074 0.8442520463
[37,] 0.2294611656 0.4589223313 0.7705388344
[38,] 0.3320227806 0.6640455613 0.6679772194
[39,] 0.5403110182 0.9193779636 0.4596889818
[40,] 0.6543874690 0.6912250620 0.3456125310
[41,] 0.7765347543 0.4469304915 0.2234652457
[42,] 0.8519916292 0.2960167416 0.1480083708
[43,] 0.8769376668 0.2461246663 0.1230623332
[44,] 0.9025796062 0.1948407876 0.0974203938
[45,] 0.9205983949 0.1588032101 0.0794016051
[46,] 0.9410732752 0.1178534497 0.0589267248
[47,] 0.9567735173 0.0864529654 0.0432264827
[48,] 0.9762057132 0.0475885736 0.0237942868
[49,] 0.9971791642 0.0056416717 0.0028208358
[50,] 0.9987779187 0.0024441625 0.0012220813
[51,] 0.9992527099 0.0014945802 0.0007472901
[52,] 0.9996895269 0.0006209463 0.0003104731
[53,] 0.9997702461 0.0004595079 0.0002297539
[54,] 0.9998110301 0.0003779399 0.0001889699
[55,] 0.9997190857 0.0005618287 0.0002809143
[56,] 0.9995955228 0.0008089543 0.0004044772
[57,] 0.9994051233 0.0011897535 0.0005948767
[58,] 0.9990383535 0.0019232929 0.0009616465
[59,] 0.9982699328 0.0034601343 0.0017300672
[60,] 0.9974744702 0.0050510597 0.0025255298
[61,] 0.9977620279 0.0044759442 0.0022379721
[62,] 0.9956457864 0.0087084273 0.0043542136
[63,] 0.9943965647 0.0112068705 0.0056034353
[64,] 0.9969715275 0.0060569450 0.0030284725
[65,] 0.9984820197 0.0030359605 0.0015179803
[66,] 0.9976990813 0.0046018374 0.0023009187
[67,] 0.9960284736 0.0079430528 0.0039715264
[68,] 0.9947096187 0.0105807626 0.0052903813
[69,] 0.9878839017 0.0242321966 0.0121160983
[70,] 0.9736936446 0.0526127108 0.0263063554
[71,] 0.9654997954 0.0690004092 0.0345002046
[72,] 0.9457571475 0.1084857049 0.0542428525
[73,] 0.8890821593 0.2218356814 0.1109178407
[74,] 0.7899753736 0.4200492528 0.2100246264
[75,] 0.6500996065 0.6998007871 0.3499003935
> postscript(file="/var/www/html/rcomp/tmp/1s1ix1229717325.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/2k7211229717325.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/353b11229717325.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/4951i1229717325.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/59vef1229717325.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 7
-3326.066 -19095.880 -36310.889 -25303.715 -44016.192 -39441.779 -3525.471
8 9 10 11 12 13 14
11022.560 -21121.289 -34075.262 -46784.799 -39462.956 -37106.725 -39507.852
15 16 17 18 19 20 21
-50820.948 -42358.302 -50481.640 -49001.858 -9312.354 7978.418 -2503.487
22 23 24 25 26 27 28
-12781.971 -18993.404 -8968.436 -10313.970 -14172.838 -24388.881 -27169.822
29 30 31 32 33 34 35
-29965.366 -28029.484 15752.650 31192.347 5532.000 -8196.088 -2340.916
36 37 38 39 40 41 42
-2342.234 1555.745 -1947.724 -11276.317 -6888.100 -6766.563 -4148.131
43 44 45 46 47 48 49
41019.979 60930.064 29900.298 14749.158 13618.715 14285.050 21352.403
50 51 52 53 54 55 56
16113.131 -2442.836 2495.402 4720.573 -7500.249 45018.746 64385.241
57 58 59 60 61 62 63
35214.299 35825.650 26186.904 27890.160 32879.883 24170.483 9925.207
64 65 66 67 68 69 70
13254.309 10726.702 2068.182 58125.359 63196.438 32636.580 32951.206
71 72 73 74 75 76 77
10988.861 12965.982 4080.834 5127.311 -20673.120 4562.134 -8531.660
78 79 80 81 82 83 84
-9319.621 33140.023 30859.869 11168.950 -14442.730 -25340.278 -19069.608
> postscript(file="/var/www/html/rcomp/tmp/6hsc01229717325.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 -3326.066 NA
1 -19095.880 -3326.066
2 -36310.889 -19095.880
3 -25303.715 -36310.889
4 -44016.192 -25303.715
5 -39441.779 -44016.192
6 -3525.471 -39441.779
7 11022.560 -3525.471
8 -21121.289 11022.560
9 -34075.262 -21121.289
10 -46784.799 -34075.262
11 -39462.956 -46784.799
12 -37106.725 -39462.956
13 -39507.852 -37106.725
14 -50820.948 -39507.852
15 -42358.302 -50820.948
16 -50481.640 -42358.302
17 -49001.858 -50481.640
18 -9312.354 -49001.858
19 7978.418 -9312.354
20 -2503.487 7978.418
21 -12781.971 -2503.487
22 -18993.404 -12781.971
23 -8968.436 -18993.404
24 -10313.970 -8968.436
25 -14172.838 -10313.970
26 -24388.881 -14172.838
27 -27169.822 -24388.881
28 -29965.366 -27169.822
29 -28029.484 -29965.366
30 15752.650 -28029.484
31 31192.347 15752.650
32 5532.000 31192.347
33 -8196.088 5532.000
34 -2340.916 -8196.088
35 -2342.234 -2340.916
36 1555.745 -2342.234
37 -1947.724 1555.745
38 -11276.317 -1947.724
39 -6888.100 -11276.317
40 -6766.563 -6888.100
41 -4148.131 -6766.563
42 41019.979 -4148.131
43 60930.064 41019.979
44 29900.298 60930.064
45 14749.158 29900.298
46 13618.715 14749.158
47 14285.050 13618.715
48 21352.403 14285.050
49 16113.131 21352.403
50 -2442.836 16113.131
51 2495.402 -2442.836
52 4720.573 2495.402
53 -7500.249 4720.573
54 45018.746 -7500.249
55 64385.241 45018.746
56 35214.299 64385.241
57 35825.650 35214.299
58 26186.904 35825.650
59 27890.160 26186.904
60 32879.883 27890.160
61 24170.483 32879.883
62 9925.207 24170.483
63 13254.309 9925.207
64 10726.702 13254.309
65 2068.182 10726.702
66 58125.359 2068.182
67 63196.438 58125.359
68 32636.580 63196.438
69 32951.206 32636.580
70 10988.861 32951.206
71 12965.982 10988.861
72 4080.834 12965.982
73 5127.311 4080.834
74 -20673.120 5127.311
75 4562.134 -20673.120
76 -8531.660 4562.134
77 -9319.621 -8531.660
78 33140.023 -9319.621
79 30859.869 33140.023
80 11168.950 30859.869
81 -14442.730 11168.950
82 -25340.278 -14442.730
83 -19069.608 -25340.278
84 NA -19069.608
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -19095.880 -3326.066
[2,] -36310.889 -19095.880
[3,] -25303.715 -36310.889
[4,] -44016.192 -25303.715
[5,] -39441.779 -44016.192
[6,] -3525.471 -39441.779
[7,] 11022.560 -3525.471
[8,] -21121.289 11022.560
[9,] -34075.262 -21121.289
[10,] -46784.799 -34075.262
[11,] -39462.956 -46784.799
[12,] -37106.725 -39462.956
[13,] -39507.852 -37106.725
[14,] -50820.948 -39507.852
[15,] -42358.302 -50820.948
[16,] -50481.640 -42358.302
[17,] -49001.858 -50481.640
[18,] -9312.354 -49001.858
[19,] 7978.418 -9312.354
[20,] -2503.487 7978.418
[21,] -12781.971 -2503.487
[22,] -18993.404 -12781.971
[23,] -8968.436 -18993.404
[24,] -10313.970 -8968.436
[25,] -14172.838 -10313.970
[26,] -24388.881 -14172.838
[27,] -27169.822 -24388.881
[28,] -29965.366 -27169.822
[29,] -28029.484 -29965.366
[30,] 15752.650 -28029.484
[31,] 31192.347 15752.650
[32,] 5532.000 31192.347
[33,] -8196.088 5532.000
[34,] -2340.916 -8196.088
[35,] -2342.234 -2340.916
[36,] 1555.745 -2342.234
[37,] -1947.724 1555.745
[38,] -11276.317 -1947.724
[39,] -6888.100 -11276.317
[40,] -6766.563 -6888.100
[41,] -4148.131 -6766.563
[42,] 41019.979 -4148.131
[43,] 60930.064 41019.979
[44,] 29900.298 60930.064
[45,] 14749.158 29900.298
[46,] 13618.715 14749.158
[47,] 14285.050 13618.715
[48,] 21352.403 14285.050
[49,] 16113.131 21352.403
[50,] -2442.836 16113.131
[51,] 2495.402 -2442.836
[52,] 4720.573 2495.402
[53,] -7500.249 4720.573
[54,] 45018.746 -7500.249
[55,] 64385.241 45018.746
[56,] 35214.299 64385.241
[57,] 35825.650 35214.299
[58,] 26186.904 35825.650
[59,] 27890.160 26186.904
[60,] 32879.883 27890.160
[61,] 24170.483 32879.883
[62,] 9925.207 24170.483
[63,] 13254.309 9925.207
[64,] 10726.702 13254.309
[65,] 2068.182 10726.702
[66,] 58125.359 2068.182
[67,] 63196.438 58125.359
[68,] 32636.580 63196.438
[69,] 32951.206 32636.580
[70,] 10988.861 32951.206
[71,] 12965.982 10988.861
[72,] 4080.834 12965.982
[73,] 5127.311 4080.834
[74,] -20673.120 5127.311
[75,] 4562.134 -20673.120
[76,] -8531.660 4562.134
[77,] -9319.621 -8531.660
[78,] 33140.023 -9319.621
[79,] 30859.869 33140.023
[80,] 11168.950 30859.869
[81,] -14442.730 11168.950
[82,] -25340.278 -14442.730
[83,] -19069.608 -25340.278
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -19095.880 -3326.066
2 -36310.889 -19095.880
3 -25303.715 -36310.889
4 -44016.192 -25303.715
5 -39441.779 -44016.192
6 -3525.471 -39441.779
7 11022.560 -3525.471
8 -21121.289 11022.560
9 -34075.262 -21121.289
10 -46784.799 -34075.262
11 -39462.956 -46784.799
12 -37106.725 -39462.956
13 -39507.852 -37106.725
14 -50820.948 -39507.852
15 -42358.302 -50820.948
16 -50481.640 -42358.302
17 -49001.858 -50481.640
18 -9312.354 -49001.858
19 7978.418 -9312.354
20 -2503.487 7978.418
21 -12781.971 -2503.487
22 -18993.404 -12781.971
23 -8968.436 -18993.404
24 -10313.970 -8968.436
25 -14172.838 -10313.970
26 -24388.881 -14172.838
27 -27169.822 -24388.881
28 -29965.366 -27169.822
29 -28029.484 -29965.366
30 15752.650 -28029.484
31 31192.347 15752.650
32 5532.000 31192.347
33 -8196.088 5532.000
34 -2340.916 -8196.088
35 -2342.234 -2340.916
36 1555.745 -2342.234
37 -1947.724 1555.745
38 -11276.317 -1947.724
39 -6888.100 -11276.317
40 -6766.563 -6888.100
41 -4148.131 -6766.563
42 41019.979 -4148.131
43 60930.064 41019.979
44 29900.298 60930.064
45 14749.158 29900.298
46 13618.715 14749.158
47 14285.050 13618.715
48 21352.403 14285.050
49 16113.131 21352.403
50 -2442.836 16113.131
51 2495.402 -2442.836
52 4720.573 2495.402
53 -7500.249 4720.573
54 45018.746 -7500.249
55 64385.241 45018.746
56 35214.299 64385.241
57 35825.650 35214.299
58 26186.904 35825.650
59 27890.160 26186.904
60 32879.883 27890.160
61 24170.483 32879.883
62 9925.207 24170.483
63 13254.309 9925.207
64 10726.702 13254.309
65 2068.182 10726.702
66 58125.359 2068.182
67 63196.438 58125.359
68 32636.580 63196.438
69 32951.206 32636.580
70 10988.861 32951.206
71 12965.982 10988.861
72 4080.834 12965.982
73 5127.311 4080.834
74 -20673.120 5127.311
75 4562.134 -20673.120
76 -8531.660 4562.134
77 -9319.621 -8531.660
78 33140.023 -9319.621
79 30859.869 33140.023
80 11168.950 30859.869
81 -14442.730 11168.950
82 -25340.278 -14442.730
83 -19069.608 -25340.278
> 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/7qz0d1229717325.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/89wex1229717325.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/9bde41229717325.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/10lrqq1229717325.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/11y91g1229717325.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/12qmfj1229717325.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/13kitk1229717325.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/14mkbc1229717325.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/15hylk1229717325.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/166dn21229717325.tab")
+ }
>
> system("convert tmp/1s1ix1229717325.ps tmp/1s1ix1229717325.png")
> system("convert tmp/2k7211229717325.ps tmp/2k7211229717325.png")
> system("convert tmp/353b11229717325.ps tmp/353b11229717325.png")
> system("convert tmp/4951i1229717325.ps tmp/4951i1229717325.png")
> system("convert tmp/59vef1229717325.ps tmp/59vef1229717325.png")
> system("convert tmp/6hsc01229717325.ps tmp/6hsc01229717325.png")
> system("convert tmp/7qz0d1229717325.ps tmp/7qz0d1229717325.png")
> system("convert tmp/89wex1229717325.ps tmp/89wex1229717325.png")
> system("convert tmp/9bde41229717325.ps tmp/9bde41229717325.png")
> system("convert tmp/10lrqq1229717325.ps tmp/10lrqq1229717325.png")
>
>
> proc.time()
user system elapsed
2.886 1.665 4.511