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(1,9700,2,9081,3,9084,4,9743,5,8587,6,9731,7,9563,8,9998,9,9437,10,10038,11,9918,12,9252,1,9737,2,9035,3,9133,4,9487,5,8700,6,9627,7,8947,8,9283,9,8829,10,9947,11,9628,12,9318,1,9605,2,8640,3,9214,4,9567,5,8547,6,9185,7,9470,8,9123,9,9278,10,10170,11,9434,12,9655,1,9429,2,8739,3,9552,4,9687,5,9019,6,9672,7,9206,8,9069,9,9788,10,10312,11,10105,12,9863,1,9656,2,9295,3,9946,4,9701,5,9049,6,10190,7,9706,8,9765,9,9893,10,9994,11,10433,12,10073,1,10112,2,9266,3,9820,4,10097,5,9115,6,10411,7,9678,8,10408,9,10153,10,10368,11,10581,12,10597,1,10680,2,9738,3,9556),dim=c(2,75),dimnames=list(c('Month','Monthly_Births'),1:75))
> y <- array(NA,dim=c(2,75),dimnames=list(c('Month','Monthly_Births'),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 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '2'
> par3 <- 'No Linear Trend'
> par2 <- 'Do not include Seasonal Dummies'
> par1 <- '2'
> #'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, 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
Monthly_Births Month
1 9700 1
2 9081 2
3 9084 3
4 9743 4
5 8587 5
6 9731 6
7 9563 7
8 9998 8
9 9437 9
10 10038 10
11 9918 11
12 9252 12
13 9737 1
14 9035 2
15 9133 3
16 9487 4
17 8700 5
18 9627 6
19 8947 7
20 9283 8
21 8829 9
22 9947 10
23 9628 11
24 9318 12
25 9605 1
26 8640 2
27 9214 3
28 9567 4
29 8547 5
30 9185 6
31 9470 7
32 9123 8
33 9278 9
34 10170 10
35 9434 11
36 9655 12
37 9429 1
38 8739 2
39 9552 3
40 9687 4
41 9019 5
42 9672 6
43 9206 7
44 9069 8
45 9788 9
46 10312 10
47 10105 11
48 9863 12
49 9656 1
50 9295 2
51 9946 3
52 9701 4
53 9049 5
54 10190 6
55 9706 7
56 9765 8
57 9893 9
58 9994 10
59 10433 11
60 10073 12
61 10112 1
62 9266 2
63 9820 3
64 10097 4
65 9115 5
66 10411 6
67 9678 7
68 10408 8
69 10153 9
70 10368 10
71 10581 11
72 10597 12
73 10680 1
74 9738 2
75 9556 3
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Month
9319.51 45.25
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-998.77 -379.15 66.48 303.85 1315.23
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9319.51 114.52 81.378 < 2e-16 ***
Month 45.25 15.85 2.855 0.00561 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 480.4 on 73 degrees of freedom
Multiple R-squared: 0.1004, Adjusted R-squared: 0.08809
F-statistic: 8.148 on 1 and 73 DF, p-value: 0.005607
> 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.6100027 0.7799947 0.3899973
[2,] 0.6933130 0.6133741 0.3066870
[3,] 0.5828534 0.8342932 0.4171466
[4,] 0.5525778 0.8948445 0.4474222
[5,] 0.4569219 0.9138439 0.5430781
[6,] 0.3893295 0.7786590 0.6106705
[7,] 0.2897189 0.5794378 0.7102811
[8,] 0.3431402 0.6862804 0.6568598
[9,] 0.3166298 0.6332597 0.6833702
[10,] 0.2811169 0.5622338 0.7188831
[11,] 0.2300761 0.4601522 0.7699239
[12,] 0.1682138 0.3364277 0.8317862
[13,] 0.2783097 0.5566194 0.7216903
[14,] 0.2188887 0.4377775 0.7811113
[15,] 0.2542234 0.5084468 0.7457766
[16,] 0.2138989 0.4277978 0.7861011
[17,] 0.3144288 0.6288577 0.6855712
[18,] 0.2950200 0.5900400 0.7049800
[19,] 0.2387437 0.4774874 0.7612563
[20,] 0.2268176 0.4536351 0.7731824
[21,] 0.1966754 0.3933508 0.8033246
[22,] 0.2770112 0.5540225 0.7229888
[23,] 0.2297086 0.4594172 0.7702914
[24,] 0.1888597 0.3777195 0.8111403
[25,] 0.3734259 0.7468517 0.6265741
[26,] 0.3502732 0.7005464 0.6497268
[27,] 0.3011603 0.6023206 0.6988397
[28,] 0.3203256 0.6406512 0.6796744
[29,] 0.3166640 0.6333281 0.6833360
[30,] 0.3512286 0.7024573 0.6487714
[31,] 0.3423251 0.6846503 0.6576749
[32,] 0.3153838 0.6307675 0.6846162
[33,] 0.2678591 0.5357182 0.7321409
[34,] 0.3622514 0.7245029 0.6377486
[35,] 0.3194486 0.6388972 0.6805514
[36,] 0.2859024 0.5718048 0.7140976
[37,] 0.3466596 0.6933192 0.6533404
[38,] 0.3072896 0.6145793 0.6927104
[39,] 0.3506181 0.7012361 0.6493819
[40,] 0.5084740 0.9830520 0.4915260
[41,] 0.4786364 0.9572728 0.5213636
[42,] 0.5200980 0.9598039 0.4799020
[43,] 0.4900086 0.9800173 0.5099914
[44,] 0.4641917 0.9283834 0.5358083
[45,] 0.4266000 0.8532000 0.5734000
[46,] 0.4076466 0.8152931 0.5923534
[47,] 0.3997171 0.7994343 0.6002829
[48,] 0.3511374 0.7022748 0.6488626
[49,] 0.5207031 0.9585938 0.4792969
[50,] 0.5259115 0.9481769 0.4740885
[51,] 0.4917646 0.9835291 0.5082354
[52,] 0.4592346 0.9184692 0.5407654
[53,] 0.4184657 0.8369313 0.5815343
[54,] 0.3757645 0.7515289 0.6242355
[55,] 0.3544090 0.7088179 0.6455910
[56,] 0.3156768 0.6313536 0.6843232
[57,] 0.3443395 0.6886790 0.6556605
[58,] 0.3592016 0.7184032 0.6407984
[59,] 0.2915045 0.5830089 0.7084955
[60,] 0.2440541 0.4881083 0.7559459
[61,] 0.5447366 0.9105267 0.4552634
[62,] 0.5114642 0.9770716 0.4885358
[63,] 0.5708639 0.8582722 0.4291361
[64,] 0.4744846 0.9489691 0.5255154
[65,] 0.3672166 0.7344333 0.6327834
[66,] 0.2426382 0.4852763 0.7573618
> postscript(file="/var/wessaorg/rcomp/tmp/17bqh1354809275.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/2hjxh1354809275.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/3740m1354809276.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/4m2jl1354809276.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/5oquy1354809276.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
335.2341087 -329.0179168 -371.2699422 242.4780323 -958.7739931 139.9739815
7 8 9 10 11 12
-73.2780440 316.4699306 -289.7820948 265.9658797 100.7138543 -610.5381712
13 14 15 16 17 18
372.2341087 -375.0179168 -322.2699422 -13.5219677 -845.7739931 35.9739815
19 20 21 22 23 24
-689.2780440 -398.5300694 -897.7820948 174.9658797 -189.2861457 -544.5381712
25 26 27 28 29 30
240.2341087 -770.0179168 -241.2699422 66.4780323 -998.7739931 -406.0260185
31 32 33 34 35 36
-166.2780440 -558.5300694 -448.7820948 397.9658797 -383.2861457 -207.5381712
37 38 39 40 41 42
64.2341087 -671.0179168 96.7300578 186.4780323 -526.7739931 80.9739815
43 44 45 46 47 48
-430.2780440 -612.5300694 61.2179052 539.9658797 287.7138543 0.4618288
49 50 51 52 53 54
291.2341087 -115.0179168 490.7300578 200.4780323 -496.7739931 598.9739815
55 56 57 58 59 60
69.7219560 83.4699306 166.2179052 221.9658797 615.7138543 210.4618288
61 62 63 64 65 66
747.2341087 -144.0179168 364.7300578 596.4780323 -430.7739931 819.9739815
67 68 69 70 71 72
41.7219560 726.4699306 426.2179052 595.9658797 763.7138543 734.4618288
73 74 75
1315.2341087 327.9820832 100.7300578
> postscript(file="/var/wessaorg/rcomp/tmp/6kvf41354809276.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 335.2341087 NA
1 -329.0179168 335.2341087
2 -371.2699422 -329.0179168
3 242.4780323 -371.2699422
4 -958.7739931 242.4780323
5 139.9739815 -958.7739931
6 -73.2780440 139.9739815
7 316.4699306 -73.2780440
8 -289.7820948 316.4699306
9 265.9658797 -289.7820948
10 100.7138543 265.9658797
11 -610.5381712 100.7138543
12 372.2341087 -610.5381712
13 -375.0179168 372.2341087
14 -322.2699422 -375.0179168
15 -13.5219677 -322.2699422
16 -845.7739931 -13.5219677
17 35.9739815 -845.7739931
18 -689.2780440 35.9739815
19 -398.5300694 -689.2780440
20 -897.7820948 -398.5300694
21 174.9658797 -897.7820948
22 -189.2861457 174.9658797
23 -544.5381712 -189.2861457
24 240.2341087 -544.5381712
25 -770.0179168 240.2341087
26 -241.2699422 -770.0179168
27 66.4780323 -241.2699422
28 -998.7739931 66.4780323
29 -406.0260185 -998.7739931
30 -166.2780440 -406.0260185
31 -558.5300694 -166.2780440
32 -448.7820948 -558.5300694
33 397.9658797 -448.7820948
34 -383.2861457 397.9658797
35 -207.5381712 -383.2861457
36 64.2341087 -207.5381712
37 -671.0179168 64.2341087
38 96.7300578 -671.0179168
39 186.4780323 96.7300578
40 -526.7739931 186.4780323
41 80.9739815 -526.7739931
42 -430.2780440 80.9739815
43 -612.5300694 -430.2780440
44 61.2179052 -612.5300694
45 539.9658797 61.2179052
46 287.7138543 539.9658797
47 0.4618288 287.7138543
48 291.2341087 0.4618288
49 -115.0179168 291.2341087
50 490.7300578 -115.0179168
51 200.4780323 490.7300578
52 -496.7739931 200.4780323
53 598.9739815 -496.7739931
54 69.7219560 598.9739815
55 83.4699306 69.7219560
56 166.2179052 83.4699306
57 221.9658797 166.2179052
58 615.7138543 221.9658797
59 210.4618288 615.7138543
60 747.2341087 210.4618288
61 -144.0179168 747.2341087
62 364.7300578 -144.0179168
63 596.4780323 364.7300578
64 -430.7739931 596.4780323
65 819.9739815 -430.7739931
66 41.7219560 819.9739815
67 726.4699306 41.7219560
68 426.2179052 726.4699306
69 595.9658797 426.2179052
70 763.7138543 595.9658797
71 734.4618288 763.7138543
72 1315.2341087 734.4618288
73 327.9820832 1315.2341087
74 100.7300578 327.9820832
75 NA 100.7300578
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -329.0179168 335.2341087
[2,] -371.2699422 -329.0179168
[3,] 242.4780323 -371.2699422
[4,] -958.7739931 242.4780323
[5,] 139.9739815 -958.7739931
[6,] -73.2780440 139.9739815
[7,] 316.4699306 -73.2780440
[8,] -289.7820948 316.4699306
[9,] 265.9658797 -289.7820948
[10,] 100.7138543 265.9658797
[11,] -610.5381712 100.7138543
[12,] 372.2341087 -610.5381712
[13,] -375.0179168 372.2341087
[14,] -322.2699422 -375.0179168
[15,] -13.5219677 -322.2699422
[16,] -845.7739931 -13.5219677
[17,] 35.9739815 -845.7739931
[18,] -689.2780440 35.9739815
[19,] -398.5300694 -689.2780440
[20,] -897.7820948 -398.5300694
[21,] 174.9658797 -897.7820948
[22,] -189.2861457 174.9658797
[23,] -544.5381712 -189.2861457
[24,] 240.2341087 -544.5381712
[25,] -770.0179168 240.2341087
[26,] -241.2699422 -770.0179168
[27,] 66.4780323 -241.2699422
[28,] -998.7739931 66.4780323
[29,] -406.0260185 -998.7739931
[30,] -166.2780440 -406.0260185
[31,] -558.5300694 -166.2780440
[32,] -448.7820948 -558.5300694
[33,] 397.9658797 -448.7820948
[34,] -383.2861457 397.9658797
[35,] -207.5381712 -383.2861457
[36,] 64.2341087 -207.5381712
[37,] -671.0179168 64.2341087
[38,] 96.7300578 -671.0179168
[39,] 186.4780323 96.7300578
[40,] -526.7739931 186.4780323
[41,] 80.9739815 -526.7739931
[42,] -430.2780440 80.9739815
[43,] -612.5300694 -430.2780440
[44,] 61.2179052 -612.5300694
[45,] 539.9658797 61.2179052
[46,] 287.7138543 539.9658797
[47,] 0.4618288 287.7138543
[48,] 291.2341087 0.4618288
[49,] -115.0179168 291.2341087
[50,] 490.7300578 -115.0179168
[51,] 200.4780323 490.7300578
[52,] -496.7739931 200.4780323
[53,] 598.9739815 -496.7739931
[54,] 69.7219560 598.9739815
[55,] 83.4699306 69.7219560
[56,] 166.2179052 83.4699306
[57,] 221.9658797 166.2179052
[58,] 615.7138543 221.9658797
[59,] 210.4618288 615.7138543
[60,] 747.2341087 210.4618288
[61,] -144.0179168 747.2341087
[62,] 364.7300578 -144.0179168
[63,] 596.4780323 364.7300578
[64,] -430.7739931 596.4780323
[65,] 819.9739815 -430.7739931
[66,] 41.7219560 819.9739815
[67,] 726.4699306 41.7219560
[68,] 426.2179052 726.4699306
[69,] 595.9658797 426.2179052
[70,] 763.7138543 595.9658797
[71,] 734.4618288 763.7138543
[72,] 1315.2341087 734.4618288
[73,] 327.9820832 1315.2341087
[74,] 100.7300578 327.9820832
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -329.0179168 335.2341087
2 -371.2699422 -329.0179168
3 242.4780323 -371.2699422
4 -958.7739931 242.4780323
5 139.9739815 -958.7739931
6 -73.2780440 139.9739815
7 316.4699306 -73.2780440
8 -289.7820948 316.4699306
9 265.9658797 -289.7820948
10 100.7138543 265.9658797
11 -610.5381712 100.7138543
12 372.2341087 -610.5381712
13 -375.0179168 372.2341087
14 -322.2699422 -375.0179168
15 -13.5219677 -322.2699422
16 -845.7739931 -13.5219677
17 35.9739815 -845.7739931
18 -689.2780440 35.9739815
19 -398.5300694 -689.2780440
20 -897.7820948 -398.5300694
21 174.9658797 -897.7820948
22 -189.2861457 174.9658797
23 -544.5381712 -189.2861457
24 240.2341087 -544.5381712
25 -770.0179168 240.2341087
26 -241.2699422 -770.0179168
27 66.4780323 -241.2699422
28 -998.7739931 66.4780323
29 -406.0260185 -998.7739931
30 -166.2780440 -406.0260185
31 -558.5300694 -166.2780440
32 -448.7820948 -558.5300694
33 397.9658797 -448.7820948
34 -383.2861457 397.9658797
35 -207.5381712 -383.2861457
36 64.2341087 -207.5381712
37 -671.0179168 64.2341087
38 96.7300578 -671.0179168
39 186.4780323 96.7300578
40 -526.7739931 186.4780323
41 80.9739815 -526.7739931
42 -430.2780440 80.9739815
43 -612.5300694 -430.2780440
44 61.2179052 -612.5300694
45 539.9658797 61.2179052
46 287.7138543 539.9658797
47 0.4618288 287.7138543
48 291.2341087 0.4618288
49 -115.0179168 291.2341087
50 490.7300578 -115.0179168
51 200.4780323 490.7300578
52 -496.7739931 200.4780323
53 598.9739815 -496.7739931
54 69.7219560 598.9739815
55 83.4699306 69.7219560
56 166.2179052 83.4699306
57 221.9658797 166.2179052
58 615.7138543 221.9658797
59 210.4618288 615.7138543
60 747.2341087 210.4618288
61 -144.0179168 747.2341087
62 364.7300578 -144.0179168
63 596.4780323 364.7300578
64 -430.7739931 596.4780323
65 819.9739815 -430.7739931
66 41.7219560 819.9739815
67 726.4699306 41.7219560
68 426.2179052 726.4699306
69 595.9658797 426.2179052
70 763.7138543 595.9658797
71 734.4618288 763.7138543
72 1315.2341087 734.4618288
73 327.9820832 1315.2341087
74 100.7300578 327.9820832
> 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/70rtk1354809276.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/8ajr01354809276.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/9ai161354809276.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/10evnd1354809276.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/11s2tq1354809276.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/12485h1354809276.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/13p4wo1354809276.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/14wp791354809276.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/155cl31354809276.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/16ekt01354809276.tab")
+ }
>
> try(system("convert tmp/17bqh1354809275.ps tmp/17bqh1354809275.png",intern=TRUE))
character(0)
> try(system("convert tmp/2hjxh1354809275.ps tmp/2hjxh1354809275.png",intern=TRUE))
character(0)
> try(system("convert tmp/3740m1354809276.ps tmp/3740m1354809276.png",intern=TRUE))
character(0)
> try(system("convert tmp/4m2jl1354809276.ps tmp/4m2jl1354809276.png",intern=TRUE))
character(0)
> try(system("convert tmp/5oquy1354809276.ps tmp/5oquy1354809276.png",intern=TRUE))
character(0)
> try(system("convert tmp/6kvf41354809276.ps tmp/6kvf41354809276.png",intern=TRUE))
character(0)
> try(system("convert tmp/70rtk1354809276.ps tmp/70rtk1354809276.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ajr01354809276.ps tmp/8ajr01354809276.png",intern=TRUE))
character(0)
> try(system("convert tmp/9ai161354809276.ps tmp/9ai161354809276.png",intern=TRUE))
character(0)
> try(system("convert tmp/10evnd1354809276.ps tmp/10evnd1354809276.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.070 1.121 7.205