R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-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(9,676,8,642,9,402,9,610,9,294,9,448,10,319,9,548,9,801,9,596,8,923,9,746,9,829,9,125,9,782,9,441,9,162,9,915,10,444,10,209,9,985,9,842,9,429,10,132,9,849,9,172,10,313,9,819,9,955,10,048,10,082,10,541,10,208,10,233,9,439,9,963,10,158,9,225,10,474,9,757,10,490,10,281,10,444,10,640,10,695,10,786,9,832,9,747,10,158,9,225,10,474,9,757,10,411,9,511,10,402,9,701,10,540,10,112,10,915,11,183,10,384,10,834,9,886,10,216,10,943,9,867,10,203,10,837,10,573,10,647,11,502,10,656,10,866,10,835,9,945,10,331),dim=c(2,76),dimnames=list(c('Monthyly','births'),1:76))
> y <- array(NA,dim=c(2,76),dimnames=list(c('Monthyly','births'),1:76))
> 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 = '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
> 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
Monthyly births t
1 9 676 1
2 8 642 2
3 9 402 3
4 9 610 4
5 9 294 5
6 9 448 6
7 10 319 7
8 9 548 8
9 9 801 9
10 9 596 10
11 8 923 11
12 9 746 12
13 9 829 13
14 9 125 14
15 9 782 15
16 9 441 16
17 9 162 17
18 9 915 18
19 10 444 19
20 10 209 20
21 9 985 21
22 9 842 22
23 9 429 23
24 10 132 24
25 9 849 25
26 9 172 26
27 10 313 27
28 9 819 28
29 9 955 29
30 10 48 30
31 10 82 31
32 10 541 32
33 10 208 33
34 10 233 34
35 9 439 35
36 9 963 36
37 10 158 37
38 9 225 38
39 10 474 39
40 9 757 40
41 10 490 41
42 10 281 42
43 10 444 43
44 10 640 44
45 10 695 45
46 10 786 46
47 9 832 47
48 9 747 48
49 10 158 49
50 9 225 50
51 10 474 51
52 9 757 52
53 10 411 53
54 9 511 54
55 10 402 55
56 9 701 56
57 10 540 57
58 10 112 58
59 10 915 59
60 11 183 60
61 10 384 61
62 10 834 62
63 9 886 63
64 10 216 64
65 10 943 65
66 9 867 66
67 10 203 67
68 10 837 68
69 10 573 69
70 10 647 70
71 11 502 71
72 10 656 72
73 10 866 73
74 10 835 74
75 9 945 75
76 10 331 76
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) births t
9.407442 -0.000941 0.016049
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.99816 -0.29768 0.08283 0.26055 0.92548
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.4074424 0.1366376 68.850 < 2e-16 ***
births -0.0009410 0.0001848 -5.093 2.66e-06 ***
t 0.0160487 0.0022956 6.991 1.08e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4369 on 73 degrees of freedom
Multiple R-squared: 0.4838, Adjusted R-squared: 0.4696
F-statistic: 34.21 on 2 and 73 DF, p-value: 3.3e-11
> 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.66346912 0.6730618 0.3365309
[2,] 0.76316462 0.4736708 0.2368354
[3,] 0.68485901 0.6302820 0.3151410
[4,] 0.56622097 0.8675581 0.4337790
[5,] 0.47900545 0.9580109 0.5209946
[6,] 0.56278008 0.8744398 0.4372199
[7,] 0.47068697 0.9413739 0.5293130
[8,] 0.39316991 0.7863398 0.6068301
[9,] 0.56195586 0.8760883 0.4380441
[10,] 0.47913228 0.9582646 0.5208677
[11,] 0.40999392 0.8199878 0.5900061
[12,] 0.39900735 0.7980147 0.6009926
[13,] 0.33950271 0.6790054 0.6604973
[14,] 0.49434758 0.9886952 0.5056524
[15,] 0.48337840 0.9667568 0.5166216
[16,] 0.40672023 0.8134405 0.5932798
[17,] 0.33649962 0.6729992 0.6635004
[18,] 0.34165881 0.6833176 0.6583412
[19,] 0.29857296 0.5971459 0.7014270
[20,] 0.24003095 0.4800619 0.7599691
[21,] 0.33883616 0.6776723 0.6611638
[22,] 0.33352588 0.6670518 0.6664741
[23,] 0.27859525 0.5571905 0.7214048
[24,] 0.22213961 0.4442792 0.7778604
[25,] 0.17553505 0.3510701 0.8244650
[26,] 0.13592972 0.2718594 0.8640703
[27,] 0.14863856 0.2972771 0.8513614
[28,] 0.11765787 0.2353157 0.8823421
[29,] 0.09270655 0.1854131 0.9072934
[30,] 0.13381508 0.2676302 0.8661849
[31,] 0.10401695 0.2080339 0.8959831
[32,] 0.07787057 0.1557411 0.9221294
[33,] 0.17001252 0.3400250 0.8299875
[34,] 0.15659286 0.3131857 0.8434071
[35,] 0.14664776 0.2932955 0.8533522
[36,] 0.13218548 0.2643710 0.8678145
[37,] 0.10241075 0.2048215 0.8975893
[38,] 0.08571651 0.1714330 0.9142835
[39,] 0.08762922 0.1752584 0.9123708
[40,] 0.09932944 0.1986589 0.9006706
[41,] 0.14202458 0.2840492 0.8579754
[42,] 0.13697524 0.2739505 0.8630248
[43,] 0.13797731 0.2759546 0.8620227
[44,] 0.10713598 0.2142720 0.8928640
[45,] 0.26385817 0.5277163 0.7361418
[46,] 0.22830331 0.4566066 0.7716967
[47,] 0.23151172 0.4630234 0.7684883
[48,] 0.18563997 0.3712799 0.8143600
[49,] 0.29323034 0.5864607 0.7067697
[50,] 0.23251153 0.4650231 0.7674885
[51,] 0.32668352 0.6533670 0.6733165
[52,] 0.26380619 0.5276124 0.7361938
[53,] 0.24847531 0.4969506 0.7515247
[54,] 0.23398707 0.4679741 0.7660129
[55,] 0.33676818 0.6735364 0.6632318
[56,] 0.26033903 0.5206781 0.7396610
[57,] 0.24872259 0.4974452 0.7512774
[58,] 0.28067532 0.5613506 0.7193247
[59,] 0.22790288 0.4558058 0.7720971
[60,] 0.20909489 0.4181898 0.7909051
[61,] 0.34270441 0.6854088 0.6572956
[62,] 0.46800216 0.9360043 0.5319978
[63,] 0.35063768 0.7012754 0.6493623
[64,] 0.38697234 0.7739447 0.6130277
[65,] 0.55876938 0.8824612 0.4412306
> postscript(file="/var/wessaorg/rcomp/tmp/1sq8b1322498792.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/wessaorg/rcomp/tmp/26oh51322498792.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/wessaorg/rcomp/tmp/3iins1322498792.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/wessaorg/rcomp/tmp/4td2o1322498792.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/wessaorg/rcomp/tmp/5cc5x1322498792.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 = 76
Frequency = 1
1 2 3 4 5 6
0.212621030 -0.835421516 -0.077308895 0.102369189 -0.211033760 -0.082169370
7 8 9 10 11 12
0.780393622 -0.020167413 0.201855417 -0.007097161 -0.715440751 0.101954513
13 14 15 16 17 18
0.164008305 -0.514500447 0.087684095 -0.249243713 -0.527829872 0.164690127
19 20 21 22 23 24
0.705433055 0.468250648 0.182413516 0.031802587 -0.372876813 0.331599131
25 26 27 28 29 30
-0.009756667 -0.662858572 0.453772891 -0.086132712 0.025793780 0.156263177
31 32 33 34 35 36
0.172208246 0.588075909 0.258676056 0.266152176 -0.556051728 -0.079019433
37 38 39 40 41 42
0.147431386 -0.805570731 0.412688121 -0.337059220 0.395646554 0.182929999
43 44 45 46 47 48
0.320263338 0.488649490 0.524355440 0.593937187 -0.378825812 -0.474859069
49 50 51 52 53 54
-0.045153472 -0.998155590 0.220103262 -0.529644078 0.128723142 -0.793226162
55 56 57 58 59 60
0.088156717 -0.646534714 0.185916459 -0.232877856 0.506691859 0.801835265
61 62 63 64 65 66
-0.025073611 0.382325103 -0.584791929 -0.231306874 0.436747271 -0.650817036
67 68 69 70 71 72
-0.291686015 0.288855657 0.024384413 0.077969256 0.925476339 0.054340729
73 74 75 76
0.235900802 0.190681239 -0.721858122 -0.315677384
> postscript(file="/var/wessaorg/rcomp/tmp/67bf91322498792.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 = 76
Frequency = 1
lag(myerror, k = 1) myerror
0 0.212621030 NA
1 -0.835421516 0.212621030
2 -0.077308895 -0.835421516
3 0.102369189 -0.077308895
4 -0.211033760 0.102369189
5 -0.082169370 -0.211033760
6 0.780393622 -0.082169370
7 -0.020167413 0.780393622
8 0.201855417 -0.020167413
9 -0.007097161 0.201855417
10 -0.715440751 -0.007097161
11 0.101954513 -0.715440751
12 0.164008305 0.101954513
13 -0.514500447 0.164008305
14 0.087684095 -0.514500447
15 -0.249243713 0.087684095
16 -0.527829872 -0.249243713
17 0.164690127 -0.527829872
18 0.705433055 0.164690127
19 0.468250648 0.705433055
20 0.182413516 0.468250648
21 0.031802587 0.182413516
22 -0.372876813 0.031802587
23 0.331599131 -0.372876813
24 -0.009756667 0.331599131
25 -0.662858572 -0.009756667
26 0.453772891 -0.662858572
27 -0.086132712 0.453772891
28 0.025793780 -0.086132712
29 0.156263177 0.025793780
30 0.172208246 0.156263177
31 0.588075909 0.172208246
32 0.258676056 0.588075909
33 0.266152176 0.258676056
34 -0.556051728 0.266152176
35 -0.079019433 -0.556051728
36 0.147431386 -0.079019433
37 -0.805570731 0.147431386
38 0.412688121 -0.805570731
39 -0.337059220 0.412688121
40 0.395646554 -0.337059220
41 0.182929999 0.395646554
42 0.320263338 0.182929999
43 0.488649490 0.320263338
44 0.524355440 0.488649490
45 0.593937187 0.524355440
46 -0.378825812 0.593937187
47 -0.474859069 -0.378825812
48 -0.045153472 -0.474859069
49 -0.998155590 -0.045153472
50 0.220103262 -0.998155590
51 -0.529644078 0.220103262
52 0.128723142 -0.529644078
53 -0.793226162 0.128723142
54 0.088156717 -0.793226162
55 -0.646534714 0.088156717
56 0.185916459 -0.646534714
57 -0.232877856 0.185916459
58 0.506691859 -0.232877856
59 0.801835265 0.506691859
60 -0.025073611 0.801835265
61 0.382325103 -0.025073611
62 -0.584791929 0.382325103
63 -0.231306874 -0.584791929
64 0.436747271 -0.231306874
65 -0.650817036 0.436747271
66 -0.291686015 -0.650817036
67 0.288855657 -0.291686015
68 0.024384413 0.288855657
69 0.077969256 0.024384413
70 0.925476339 0.077969256
71 0.054340729 0.925476339
72 0.235900802 0.054340729
73 0.190681239 0.235900802
74 -0.721858122 0.190681239
75 -0.315677384 -0.721858122
76 NA -0.315677384
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.835421516 0.212621030
[2,] -0.077308895 -0.835421516
[3,] 0.102369189 -0.077308895
[4,] -0.211033760 0.102369189
[5,] -0.082169370 -0.211033760
[6,] 0.780393622 -0.082169370
[7,] -0.020167413 0.780393622
[8,] 0.201855417 -0.020167413
[9,] -0.007097161 0.201855417
[10,] -0.715440751 -0.007097161
[11,] 0.101954513 -0.715440751
[12,] 0.164008305 0.101954513
[13,] -0.514500447 0.164008305
[14,] 0.087684095 -0.514500447
[15,] -0.249243713 0.087684095
[16,] -0.527829872 -0.249243713
[17,] 0.164690127 -0.527829872
[18,] 0.705433055 0.164690127
[19,] 0.468250648 0.705433055
[20,] 0.182413516 0.468250648
[21,] 0.031802587 0.182413516
[22,] -0.372876813 0.031802587
[23,] 0.331599131 -0.372876813
[24,] -0.009756667 0.331599131
[25,] -0.662858572 -0.009756667
[26,] 0.453772891 -0.662858572
[27,] -0.086132712 0.453772891
[28,] 0.025793780 -0.086132712
[29,] 0.156263177 0.025793780
[30,] 0.172208246 0.156263177
[31,] 0.588075909 0.172208246
[32,] 0.258676056 0.588075909
[33,] 0.266152176 0.258676056
[34,] -0.556051728 0.266152176
[35,] -0.079019433 -0.556051728
[36,] 0.147431386 -0.079019433
[37,] -0.805570731 0.147431386
[38,] 0.412688121 -0.805570731
[39,] -0.337059220 0.412688121
[40,] 0.395646554 -0.337059220
[41,] 0.182929999 0.395646554
[42,] 0.320263338 0.182929999
[43,] 0.488649490 0.320263338
[44,] 0.524355440 0.488649490
[45,] 0.593937187 0.524355440
[46,] -0.378825812 0.593937187
[47,] -0.474859069 -0.378825812
[48,] -0.045153472 -0.474859069
[49,] -0.998155590 -0.045153472
[50,] 0.220103262 -0.998155590
[51,] -0.529644078 0.220103262
[52,] 0.128723142 -0.529644078
[53,] -0.793226162 0.128723142
[54,] 0.088156717 -0.793226162
[55,] -0.646534714 0.088156717
[56,] 0.185916459 -0.646534714
[57,] -0.232877856 0.185916459
[58,] 0.506691859 -0.232877856
[59,] 0.801835265 0.506691859
[60,] -0.025073611 0.801835265
[61,] 0.382325103 -0.025073611
[62,] -0.584791929 0.382325103
[63,] -0.231306874 -0.584791929
[64,] 0.436747271 -0.231306874
[65,] -0.650817036 0.436747271
[66,] -0.291686015 -0.650817036
[67,] 0.288855657 -0.291686015
[68,] 0.024384413 0.288855657
[69,] 0.077969256 0.024384413
[70,] 0.925476339 0.077969256
[71,] 0.054340729 0.925476339
[72,] 0.235900802 0.054340729
[73,] 0.190681239 0.235900802
[74,] -0.721858122 0.190681239
[75,] -0.315677384 -0.721858122
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.835421516 0.212621030
2 -0.077308895 -0.835421516
3 0.102369189 -0.077308895
4 -0.211033760 0.102369189
5 -0.082169370 -0.211033760
6 0.780393622 -0.082169370
7 -0.020167413 0.780393622
8 0.201855417 -0.020167413
9 -0.007097161 0.201855417
10 -0.715440751 -0.007097161
11 0.101954513 -0.715440751
12 0.164008305 0.101954513
13 -0.514500447 0.164008305
14 0.087684095 -0.514500447
15 -0.249243713 0.087684095
16 -0.527829872 -0.249243713
17 0.164690127 -0.527829872
18 0.705433055 0.164690127
19 0.468250648 0.705433055
20 0.182413516 0.468250648
21 0.031802587 0.182413516
22 -0.372876813 0.031802587
23 0.331599131 -0.372876813
24 -0.009756667 0.331599131
25 -0.662858572 -0.009756667
26 0.453772891 -0.662858572
27 -0.086132712 0.453772891
28 0.025793780 -0.086132712
29 0.156263177 0.025793780
30 0.172208246 0.156263177
31 0.588075909 0.172208246
32 0.258676056 0.588075909
33 0.266152176 0.258676056
34 -0.556051728 0.266152176
35 -0.079019433 -0.556051728
36 0.147431386 -0.079019433
37 -0.805570731 0.147431386
38 0.412688121 -0.805570731
39 -0.337059220 0.412688121
40 0.395646554 -0.337059220
41 0.182929999 0.395646554
42 0.320263338 0.182929999
43 0.488649490 0.320263338
44 0.524355440 0.488649490
45 0.593937187 0.524355440
46 -0.378825812 0.593937187
47 -0.474859069 -0.378825812
48 -0.045153472 -0.474859069
49 -0.998155590 -0.045153472
50 0.220103262 -0.998155590
51 -0.529644078 0.220103262
52 0.128723142 -0.529644078
53 -0.793226162 0.128723142
54 0.088156717 -0.793226162
55 -0.646534714 0.088156717
56 0.185916459 -0.646534714
57 -0.232877856 0.185916459
58 0.506691859 -0.232877856
59 0.801835265 0.506691859
60 -0.025073611 0.801835265
61 0.382325103 -0.025073611
62 -0.584791929 0.382325103
63 -0.231306874 -0.584791929
64 0.436747271 -0.231306874
65 -0.650817036 0.436747271
66 -0.291686015 -0.650817036
67 0.288855657 -0.291686015
68 0.024384413 0.288855657
69 0.077969256 0.024384413
70 0.925476339 0.077969256
71 0.054340729 0.925476339
72 0.235900802 0.054340729
73 0.190681239 0.235900802
74 -0.721858122 0.190681239
75 -0.315677384 -0.721858122
> 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/wessaorg/rcomp/tmp/7ilyn1322498792.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/wessaorg/rcomp/tmp/8navg1322498792.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/wessaorg/rcomp/tmp/903s61322498792.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/wessaorg/rcomp/tmp/10s9ca1322498792.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11s02l1322498792.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/wessaorg/rcomp/tmp/12d59q1322498792.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/wessaorg/rcomp/tmp/13ieop1322498792.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/wessaorg/rcomp/tmp/14p4z81322498792.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/wessaorg/rcomp/tmp/15pxwd1322498792.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/wessaorg/rcomp/tmp/167ru81322498792.tab")
+ }
>
> try(system("convert tmp/1sq8b1322498792.ps tmp/1sq8b1322498792.png",intern=TRUE))
character(0)
> try(system("convert tmp/26oh51322498792.ps tmp/26oh51322498792.png",intern=TRUE))
character(0)
> try(system("convert tmp/3iins1322498792.ps tmp/3iins1322498792.png",intern=TRUE))
character(0)
> try(system("convert tmp/4td2o1322498792.ps tmp/4td2o1322498792.png",intern=TRUE))
character(0)
> try(system("convert tmp/5cc5x1322498792.ps tmp/5cc5x1322498792.png",intern=TRUE))
character(0)
> try(system("convert tmp/67bf91322498792.ps tmp/67bf91322498792.png",intern=TRUE))
character(0)
> try(system("convert tmp/7ilyn1322498792.ps tmp/7ilyn1322498792.png",intern=TRUE))
character(0)
> try(system("convert tmp/8navg1322498792.ps tmp/8navg1322498792.png",intern=TRUE))
character(0)
> try(system("convert tmp/903s61322498792.ps tmp/903s61322498792.png",intern=TRUE))
character(0)
> try(system("convert tmp/10s9ca1322498792.ps tmp/10s9ca1322498792.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.370 0.538 3.936