R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-pc-linux-gnu (32-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(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('Monthlybirths'),1:75))
> y <- array(NA,dim=c(1,75),dimnames=list(c('Monthlybirths'),1:75))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par20 = ''
> par19 = ''
> par18 = ''
> par17 = ''
> par16 = ''
> par15 = ''
> par14 = ''
> par13 = ''
> par12 = ''
> par11 = ''
> par10 = ''
> par9 = ''
> par8 = ''
> par7 = ''
> par6 = ''
> par5 = ''
> par4 = ''
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from 'package:base':
as.Date, 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
Monthlybirths 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/fisher/rcomp/tmp/16j9c1355074329.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/fisher/rcomp/tmp/26p5s1355074329.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/fisher/rcomp/tmp/39eva1355074329.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/fisher/rcomp/tmp/4gxxr1355074329.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/fisher/rcomp/tmp/5yhkf1355074329.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/fisher/rcomp/tmp/6yfiu1355074329.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/fisher/rcomp/tmp/7bffs1355074329.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/fisher/rcomp/tmp/8huvs1355074329.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/fisher/rcomp/tmp/9fpzm1355074329.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/fisher/rcomp/tmp/10a2om1355074329.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/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/fisher/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/fisher/rcomp/tmp/11zt4s1355074329.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/fisher/rcomp/tmp/12cu0p1355074329.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/fisher/rcomp/tmp/13ws0f1355074329.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/fisher/rcomp/tmp/146b7n1355074330.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/fisher/rcomp/tmp/1516on1355074330.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/fisher/rcomp/tmp/161p3n1355074330.tab")
+ }
>
> try(system("convert tmp/16j9c1355074329.ps tmp/16j9c1355074329.png",intern=TRUE))
character(0)
> try(system("convert tmp/26p5s1355074329.ps tmp/26p5s1355074329.png",intern=TRUE))
character(0)
> try(system("convert tmp/39eva1355074329.ps tmp/39eva1355074329.png",intern=TRUE))
character(0)
> try(system("convert tmp/4gxxr1355074329.ps tmp/4gxxr1355074329.png",intern=TRUE))
character(0)
> try(system("convert tmp/5yhkf1355074329.ps tmp/5yhkf1355074329.png",intern=TRUE))
character(0)
> try(system("convert tmp/6yfiu1355074329.ps tmp/6yfiu1355074329.png",intern=TRUE))
character(0)
> try(system("convert tmp/7bffs1355074329.ps tmp/7bffs1355074329.png",intern=TRUE))
character(0)
> try(system("convert tmp/8huvs1355074329.ps tmp/8huvs1355074329.png",intern=TRUE))
character(0)
> try(system("convert tmp/9fpzm1355074329.ps tmp/9fpzm1355074329.png",intern=TRUE))
character(0)
> try(system("convert tmp/10a2om1355074329.ps tmp/10a2om1355074329.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.276 1.648 7.926