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.
Natural language support but running in an English locale
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(9700,9081,9084,9743,8587,9731,9563,9998,9437,10038,9918,9252,9737,9035,9133,9487,8700,9627,8947,9283,8829,9947,9628,9318,9605,8640,9214,9567,8547,9185,9470,9123,9278,10170,9434,9655,9429,8739,9552,9687,9019,9672,9206,9069,9788,10312,10105,9863,9656,9295,9946,9701,9049,10190,9706,9765,9893,9994,10433,10073,10112,9266,9820,10097,9115,10411,9678,10408,10153,10368,10581,10597,10680,9738,9556),dim=c(1,75),dimnames=list(c('Geboortes_per_maand'),1:75))
> y <- array(NA,dim=c(1,75),dimnames=list(c('Geboortes_per_maand'),1:75))
> 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
Geboortes_per_maand M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 9700 1 0 0 0 0 0 0 0 0 0 0 1
2 9081 0 1 0 0 0 0 0 0 0 0 0 2
3 9084 0 0 1 0 0 0 0 0 0 0 0 3
4 9743 0 0 0 1 0 0 0 0 0 0 0 4
5 8587 0 0 0 0 1 0 0 0 0 0 0 5
6 9731 0 0 0 0 0 1 0 0 0 0 0 6
7 9563 0 0 0 0 0 0 1 0 0 0 0 7
8 9998 0 0 0 0 0 0 0 1 0 0 0 8
9 9437 0 0 0 0 0 0 0 0 1 0 0 9
10 10038 0 0 0 0 0 0 0 0 0 1 0 10
11 9918 0 0 0 0 0 0 0 0 0 0 1 11
12 9252 0 0 0 0 0 0 0 0 0 0 0 12
13 9737 1 0 0 0 0 0 0 0 0 0 0 13
14 9035 0 1 0 0 0 0 0 0 0 0 0 14
15 9133 0 0 1 0 0 0 0 0 0 0 0 15
16 9487 0 0 0 1 0 0 0 0 0 0 0 16
17 8700 0 0 0 0 1 0 0 0 0 0 0 17
18 9627 0 0 0 0 0 1 0 0 0 0 0 18
19 8947 0 0 0 0 0 0 1 0 0 0 0 19
20 9283 0 0 0 0 0 0 0 1 0 0 0 20
21 8829 0 0 0 0 0 0 0 0 1 0 0 21
22 9947 0 0 0 0 0 0 0 0 0 1 0 22
23 9628 0 0 0 0 0 0 0 0 0 0 1 23
24 9318 0 0 0 0 0 0 0 0 0 0 0 24
25 9605 1 0 0 0 0 0 0 0 0 0 0 25
26 8640 0 1 0 0 0 0 0 0 0 0 0 26
27 9214 0 0 1 0 0 0 0 0 0 0 0 27
28 9567 0 0 0 1 0 0 0 0 0 0 0 28
29 8547 0 0 0 0 1 0 0 0 0 0 0 29
30 9185 0 0 0 0 0 1 0 0 0 0 0 30
31 9470 0 0 0 0 0 0 1 0 0 0 0 31
32 9123 0 0 0 0 0 0 0 1 0 0 0 32
33 9278 0 0 0 0 0 0 0 0 1 0 0 33
34 10170 0 0 0 0 0 0 0 0 0 1 0 34
35 9434 0 0 0 0 0 0 0 0 0 0 1 35
36 9655 0 0 0 0 0 0 0 0 0 0 0 36
37 9429 1 0 0 0 0 0 0 0 0 0 0 37
38 8739 0 1 0 0 0 0 0 0 0 0 0 38
39 9552 0 0 1 0 0 0 0 0 0 0 0 39
40 9687 0 0 0 1 0 0 0 0 0 0 0 40
41 9019 0 0 0 0 1 0 0 0 0 0 0 41
42 9672 0 0 0 0 0 1 0 0 0 0 0 42
43 9206 0 0 0 0 0 0 1 0 0 0 0 43
44 9069 0 0 0 0 0 0 0 1 0 0 0 44
45 9788 0 0 0 0 0 0 0 0 1 0 0 45
46 10312 0 0 0 0 0 0 0 0 0 1 0 46
47 10105 0 0 0 0 0 0 0 0 0 0 1 47
48 9863 0 0 0 0 0 0 0 0 0 0 0 48
49 9656 1 0 0 0 0 0 0 0 0 0 0 49
50 9295 0 1 0 0 0 0 0 0 0 0 0 50
51 9946 0 0 1 0 0 0 0 0 0 0 0 51
52 9701 0 0 0 1 0 0 0 0 0 0 0 52
53 9049 0 0 0 0 1 0 0 0 0 0 0 53
54 10190 0 0 0 0 0 1 0 0 0 0 0 54
55 9706 0 0 0 0 0 0 1 0 0 0 0 55
56 9765 0 0 0 0 0 0 0 1 0 0 0 56
57 9893 0 0 0 0 0 0 0 0 1 0 0 57
58 9994 0 0 0 0 0 0 0 0 0 1 0 58
59 10433 0 0 0 0 0 0 0 0 0 0 1 59
60 10073 0 0 0 0 0 0 0 0 0 0 0 60
61 10112 1 0 0 0 0 0 0 0 0 0 0 61
62 9266 0 1 0 0 0 0 0 0 0 0 0 62
63 9820 0 0 1 0 0 0 0 0 0 0 0 63
64 10097 0 0 0 1 0 0 0 0 0 0 0 64
65 9115 0 0 0 0 1 0 0 0 0 0 0 65
66 10411 0 0 0 0 0 1 0 0 0 0 0 66
67 9678 0 0 0 0 0 0 1 0 0 0 0 67
68 10408 0 0 0 0 0 0 0 1 0 0 0 68
69 10153 0 0 0 0 0 0 0 0 1 0 0 69
70 10368 0 0 0 0 0 0 0 0 0 1 0 70
71 10581 0 0 0 0 0 0 0 0 0 0 1 71
72 10597 0 0 0 0 0 0 0 0 0 0 0 72
73 10680 1 0 0 0 0 0 0 0 0 0 0 73
74 9738 0 1 0 0 0 0 0 0 0 0 0 74
75 9556 0 0 1 0 0 0 0 0 0 0 0 75
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
9330.587 107.621 -635.532 -287.828 8.745 -879.764
M6 M7 M8 M9 M10 M11
75.726 -309.617 -141.294 -196.970 367.186 234.510
t
11.010
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-604.73 -193.52 14.66 187.48 720.63
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9330.587 136.432 68.390 < 2e-16 ***
M1 107.621 162.985 0.660 0.511500
M2 -635.532 162.917 -3.901 0.000238 ***
M3 -287.828 162.864 -1.767 0.082101 .
M4 8.745 169.407 0.052 0.958995
M5 -879.764 169.298 -5.197 2.40e-06 ***
M6 75.726 169.203 0.448 0.656043
M7 -309.617 169.123 -1.831 0.071949 .
M8 -141.294 169.058 -0.836 0.406492
M9 -196.970 169.007 -1.165 0.248298
M10 367.186 168.970 2.173 0.033604 *
M11 234.510 168.949 1.388 0.170088
t 11.010 1.569 7.017 2.01e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 292.6 on 62 degrees of freedom
Multiple R-squared: 0.7165, Adjusted R-squared: 0.6617
F-statistic: 13.06 on 12 and 62 DF, p-value: 8.078e-13
> 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.08579503 0.17159006 0.91420497
[2,] 0.05249914 0.10499828 0.94750086
[3,] 0.02310440 0.04620879 0.97689560
[4,] 0.23149521 0.46299042 0.76850479
[5,] 0.45579254 0.91158507 0.54420746
[6,] 0.50462485 0.99075029 0.49537515
[7,] 0.43767767 0.87535533 0.56232233
[8,] 0.34035187 0.68070374 0.65964813
[9,] 0.31085586 0.62171172 0.68914414
[10,] 0.26723634 0.53447267 0.73276366
[11,] 0.20684988 0.41369977 0.79315012
[12,] 0.23991625 0.47983251 0.76008375
[13,] 0.20563145 0.41126291 0.79436855
[14,] 0.15326896 0.30653792 0.84673104
[15,] 0.17300846 0.34601692 0.82699154
[16,] 0.26800300 0.53600600 0.73199700
[17,] 0.26141642 0.52283284 0.73858358
[18,] 0.26845892 0.53691783 0.73154108
[19,] 0.34299513 0.68599026 0.65700487
[20,] 0.35877367 0.71754735 0.64122633
[21,] 0.43190572 0.86381144 0.56809428
[22,] 0.38036295 0.76072589 0.61963705
[23,] 0.32910088 0.65820175 0.67089912
[24,] 0.45050500 0.90101001 0.54949500
[25,] 0.40319020 0.80638040 0.59680980
[26,] 0.48509413 0.97018826 0.51490587
[27,] 0.45431710 0.90863420 0.54568290
[28,] 0.38075627 0.76151254 0.61924373
[29,] 0.60612685 0.78774629 0.39387315
[30,] 0.67499677 0.65000646 0.32500323
[31,] 0.75229077 0.49541846 0.24770923
[32,] 0.72881192 0.54237616 0.27118808
[33,] 0.70413873 0.59172254 0.29586127
[34,] 0.74712683 0.50574635 0.25287317
[35,] 0.69942199 0.60115602 0.30057801
[36,] 0.91622456 0.16755088 0.08377544
[37,] 0.87268706 0.25462588 0.12731294
[38,] 0.83758669 0.32482661 0.16241331
[39,] 0.79132071 0.41735857 0.20867929
[40,] 0.78948386 0.42103228 0.21051614
[41,] 0.77592534 0.44814933 0.22407466
[42,] 0.67305864 0.65388272 0.32694136
[43,] 0.53945231 0.92109538 0.46054769
[44,] 0.41611044 0.83222088 0.58388956
> postscript(file="/var/www/html/freestat/rcomp/tmp/1spft1291114326.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/html/freestat/rcomp/tmp/22yee1291114326.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/html/freestat/rcomp/tmp/32yee1291114326.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/html/freestat/rcomp/tmp/42yee1291114326.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/html/freestat/rcomp/tmp/52yee1291114326.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 = 75
Frequency = 1
1 2 3 4 5 6
250.782609 363.925466 8.211180 359.628364 81.128364 258.628364
7 8 9 10 11 12
464.961698 720.628364 204.295031 230.128364 231.795031 -210.704969
13 14 15 16 17 18
155.664596 185.807453 -74.906832 -28.489648 62.010352 22.510352
19 20 21 22 23 24
-283.156315 -126.489648 -535.822981 7.010352 -190.322981 -276.822981
25 26 27 28 29 30
-108.453416 -341.310559 -126.024845 -80.607660 -223.107660 -551.607660
31 32 33 34 35 36
107.725673 -418.607660 -218.940994 97.892340 -516.440994 -71.940994
37 38 39 40 41 42
-416.571429 -374.428571 79.857143 -92.725673 116.774327 -196.725673
43 44 45 46 47 48
-288.392340 -604.725673 158.940994 107.774327 22.440994 3.940994
49 50 51 52 53 54
-321.689441 49.453416 341.739130 -210.843685 14.656315 189.156315
55 56 57 58 59 60
79.489648 -40.843685 131.822981 -342.343685 218.322981 81.822981
61 62 63 64 65 66
2.192547 -111.664596 83.621118 53.038302 -51.461698 278.038302
67 68 69 70 71 72
-80.628364 470.038302 259.704969 -100.461698 234.204969 473.704969
73 74 75
438.074534 228.217391 -312.496894
> postscript(file="/var/www/html/freestat/rcomp/tmp/6v8vh1291114326.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 = 75
Frequency = 1
lag(myerror, k = 1) myerror
0 250.782609 NA
1 363.925466 250.782609
2 8.211180 363.925466
3 359.628364 8.211180
4 81.128364 359.628364
5 258.628364 81.128364
6 464.961698 258.628364
7 720.628364 464.961698
8 204.295031 720.628364
9 230.128364 204.295031
10 231.795031 230.128364
11 -210.704969 231.795031
12 155.664596 -210.704969
13 185.807453 155.664596
14 -74.906832 185.807453
15 -28.489648 -74.906832
16 62.010352 -28.489648
17 22.510352 62.010352
18 -283.156315 22.510352
19 -126.489648 -283.156315
20 -535.822981 -126.489648
21 7.010352 -535.822981
22 -190.322981 7.010352
23 -276.822981 -190.322981
24 -108.453416 -276.822981
25 -341.310559 -108.453416
26 -126.024845 -341.310559
27 -80.607660 -126.024845
28 -223.107660 -80.607660
29 -551.607660 -223.107660
30 107.725673 -551.607660
31 -418.607660 107.725673
32 -218.940994 -418.607660
33 97.892340 -218.940994
34 -516.440994 97.892340
35 -71.940994 -516.440994
36 -416.571429 -71.940994
37 -374.428571 -416.571429
38 79.857143 -374.428571
39 -92.725673 79.857143
40 116.774327 -92.725673
41 -196.725673 116.774327
42 -288.392340 -196.725673
43 -604.725673 -288.392340
44 158.940994 -604.725673
45 107.774327 158.940994
46 22.440994 107.774327
47 3.940994 22.440994
48 -321.689441 3.940994
49 49.453416 -321.689441
50 341.739130 49.453416
51 -210.843685 341.739130
52 14.656315 -210.843685
53 189.156315 14.656315
54 79.489648 189.156315
55 -40.843685 79.489648
56 131.822981 -40.843685
57 -342.343685 131.822981
58 218.322981 -342.343685
59 81.822981 218.322981
60 2.192547 81.822981
61 -111.664596 2.192547
62 83.621118 -111.664596
63 53.038302 83.621118
64 -51.461698 53.038302
65 278.038302 -51.461698
66 -80.628364 278.038302
67 470.038302 -80.628364
68 259.704969 470.038302
69 -100.461698 259.704969
70 234.204969 -100.461698
71 473.704969 234.204969
72 438.074534 473.704969
73 228.217391 438.074534
74 -312.496894 228.217391
75 NA -312.496894
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 363.925466 250.782609
[2,] 8.211180 363.925466
[3,] 359.628364 8.211180
[4,] 81.128364 359.628364
[5,] 258.628364 81.128364
[6,] 464.961698 258.628364
[7,] 720.628364 464.961698
[8,] 204.295031 720.628364
[9,] 230.128364 204.295031
[10,] 231.795031 230.128364
[11,] -210.704969 231.795031
[12,] 155.664596 -210.704969
[13,] 185.807453 155.664596
[14,] -74.906832 185.807453
[15,] -28.489648 -74.906832
[16,] 62.010352 -28.489648
[17,] 22.510352 62.010352
[18,] -283.156315 22.510352
[19,] -126.489648 -283.156315
[20,] -535.822981 -126.489648
[21,] 7.010352 -535.822981
[22,] -190.322981 7.010352
[23,] -276.822981 -190.322981
[24,] -108.453416 -276.822981
[25,] -341.310559 -108.453416
[26,] -126.024845 -341.310559
[27,] -80.607660 -126.024845
[28,] -223.107660 -80.607660
[29,] -551.607660 -223.107660
[30,] 107.725673 -551.607660
[31,] -418.607660 107.725673
[32,] -218.940994 -418.607660
[33,] 97.892340 -218.940994
[34,] -516.440994 97.892340
[35,] -71.940994 -516.440994
[36,] -416.571429 -71.940994
[37,] -374.428571 -416.571429
[38,] 79.857143 -374.428571
[39,] -92.725673 79.857143
[40,] 116.774327 -92.725673
[41,] -196.725673 116.774327
[42,] -288.392340 -196.725673
[43,] -604.725673 -288.392340
[44,] 158.940994 -604.725673
[45,] 107.774327 158.940994
[46,] 22.440994 107.774327
[47,] 3.940994 22.440994
[48,] -321.689441 3.940994
[49,] 49.453416 -321.689441
[50,] 341.739130 49.453416
[51,] -210.843685 341.739130
[52,] 14.656315 -210.843685
[53,] 189.156315 14.656315
[54,] 79.489648 189.156315
[55,] -40.843685 79.489648
[56,] 131.822981 -40.843685
[57,] -342.343685 131.822981
[58,] 218.322981 -342.343685
[59,] 81.822981 218.322981
[60,] 2.192547 81.822981
[61,] -111.664596 2.192547
[62,] 83.621118 -111.664596
[63,] 53.038302 83.621118
[64,] -51.461698 53.038302
[65,] 278.038302 -51.461698
[66,] -80.628364 278.038302
[67,] 470.038302 -80.628364
[68,] 259.704969 470.038302
[69,] -100.461698 259.704969
[70,] 234.204969 -100.461698
[71,] 473.704969 234.204969
[72,] 438.074534 473.704969
[73,] 228.217391 438.074534
[74,] -312.496894 228.217391
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 363.925466 250.782609
2 8.211180 363.925466
3 359.628364 8.211180
4 81.128364 359.628364
5 258.628364 81.128364
6 464.961698 258.628364
7 720.628364 464.961698
8 204.295031 720.628364
9 230.128364 204.295031
10 231.795031 230.128364
11 -210.704969 231.795031
12 155.664596 -210.704969
13 185.807453 155.664596
14 -74.906832 185.807453
15 -28.489648 -74.906832
16 62.010352 -28.489648
17 22.510352 62.010352
18 -283.156315 22.510352
19 -126.489648 -283.156315
20 -535.822981 -126.489648
21 7.010352 -535.822981
22 -190.322981 7.010352
23 -276.822981 -190.322981
24 -108.453416 -276.822981
25 -341.310559 -108.453416
26 -126.024845 -341.310559
27 -80.607660 -126.024845
28 -223.107660 -80.607660
29 -551.607660 -223.107660
30 107.725673 -551.607660
31 -418.607660 107.725673
32 -218.940994 -418.607660
33 97.892340 -218.940994
34 -516.440994 97.892340
35 -71.940994 -516.440994
36 -416.571429 -71.940994
37 -374.428571 -416.571429
38 79.857143 -374.428571
39 -92.725673 79.857143
40 116.774327 -92.725673
41 -196.725673 116.774327
42 -288.392340 -196.725673
43 -604.725673 -288.392340
44 158.940994 -604.725673
45 107.774327 158.940994
46 22.440994 107.774327
47 3.940994 22.440994
48 -321.689441 3.940994
49 49.453416 -321.689441
50 341.739130 49.453416
51 -210.843685 341.739130
52 14.656315 -210.843685
53 189.156315 14.656315
54 79.489648 189.156315
55 -40.843685 79.489648
56 131.822981 -40.843685
57 -342.343685 131.822981
58 218.322981 -342.343685
59 81.822981 218.322981
60 2.192547 81.822981
61 -111.664596 2.192547
62 83.621118 -111.664596
63 53.038302 83.621118
64 -51.461698 53.038302
65 278.038302 -51.461698
66 -80.628364 278.038302
67 470.038302 -80.628364
68 259.704969 470.038302
69 -100.461698 259.704969
70 234.204969 -100.461698
71 473.704969 234.204969
72 438.074534 473.704969
73 228.217391 438.074534
74 -312.496894 228.217391
> 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/7ohdk1291114326.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/html/freestat/rcomp/tmp/8ohdk1291114326.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/html/freestat/rcomp/tmp/9ohdk1291114326.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/html/freestat/rcomp/tmp/10y8u51291114326.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/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/11k9ss1291114326.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/12n9rg1291114326.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/13jj771291114326.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/144jnv1291114326.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/15qk3j1291114326.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/16t2kp1291114326.tab")
+ }
> try(system("convert tmp/1spft1291114326.ps tmp/1spft1291114326.png",intern=TRUE))
character(0)
> try(system("convert tmp/22yee1291114326.ps tmp/22yee1291114326.png",intern=TRUE))
character(0)
> try(system("convert tmp/32yee1291114326.ps tmp/32yee1291114326.png",intern=TRUE))
character(0)
> try(system("convert tmp/42yee1291114326.ps tmp/42yee1291114326.png",intern=TRUE))
character(0)
> try(system("convert tmp/52yee1291114326.ps tmp/52yee1291114326.png",intern=TRUE))
character(0)
> try(system("convert tmp/6v8vh1291114326.ps tmp/6v8vh1291114326.png",intern=TRUE))
character(0)
> try(system("convert tmp/7ohdk1291114326.ps tmp/7ohdk1291114326.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ohdk1291114326.ps tmp/8ohdk1291114326.png",intern=TRUE))
character(0)
> try(system("convert tmp/9ohdk1291114326.ps tmp/9ohdk1291114326.png",intern=TRUE))
character(0)
> try(system("convert tmp/10y8u51291114326.ps tmp/10y8u51291114326.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.035 2.517 4.360