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('births'),1:75))
> y <- array(NA,dim=c(1,75),dimnames=list(c('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 = 'Include Monthly Dummies'
> par1 = '1'
> par3 <- 'No 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, 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
births M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 9700 1 0 0 0 0 0 0 0 0 0 0
2 9081 0 1 0 0 0 0 0 0 0 0 0
3 9084 0 0 1 0 0 0 0 0 0 0 0
4 9743 0 0 0 1 0 0 0 0 0 0 0
5 8587 0 0 0 0 1 0 0 0 0 0 0
6 9731 0 0 0 0 0 1 0 0 0 0 0
7 9563 0 0 0 0 0 0 1 0 0 0 0
8 9998 0 0 0 0 0 0 0 1 0 0 0
9 9437 0 0 0 0 0 0 0 0 1 0 0
10 10038 0 0 0 0 0 0 0 0 0 1 0
11 9918 0 0 0 0 0 0 0 0 0 0 1
12 9252 0 0 0 0 0 0 0 0 0 0 0
13 9737 1 0 0 0 0 0 0 0 0 0 0
14 9035 0 1 0 0 0 0 0 0 0 0 0
15 9133 0 0 1 0 0 0 0 0 0 0 0
16 9487 0 0 0 1 0 0 0 0 0 0 0
17 8700 0 0 0 0 1 0 0 0 0 0 0
18 9627 0 0 0 0 0 1 0 0 0 0 0
19 8947 0 0 0 0 0 0 1 0 0 0 0
20 9283 0 0 0 0 0 0 0 1 0 0 0
21 8829 0 0 0 0 0 0 0 0 1 0 0
22 9947 0 0 0 0 0 0 0 0 0 1 0
23 9628 0 0 0 0 0 0 0 0 0 0 1
24 9318 0 0 0 0 0 0 0 0 0 0 0
25 9605 1 0 0 0 0 0 0 0 0 0 0
26 8640 0 1 0 0 0 0 0 0 0 0 0
27 9214 0 0 1 0 0 0 0 0 0 0 0
28 9567 0 0 0 1 0 0 0 0 0 0 0
29 8547 0 0 0 0 1 0 0 0 0 0 0
30 9185 0 0 0 0 0 1 0 0 0 0 0
31 9470 0 0 0 0 0 0 1 0 0 0 0
32 9123 0 0 0 0 0 0 0 1 0 0 0
33 9278 0 0 0 0 0 0 0 0 1 0 0
34 10170 0 0 0 0 0 0 0 0 0 1 0
35 9434 0 0 0 0 0 0 0 0 0 0 1
36 9655 0 0 0 0 0 0 0 0 0 0 0
37 9429 1 0 0 0 0 0 0 0 0 0 0
38 8739 0 1 0 0 0 0 0 0 0 0 0
39 9552 0 0 1 0 0 0 0 0 0 0 0
40 9687 0 0 0 1 0 0 0 0 0 0 0
41 9019 0 0 0 0 1 0 0 0 0 0 0
42 9672 0 0 0 0 0 1 0 0 0 0 0
43 9206 0 0 0 0 0 0 1 0 0 0 0
44 9069 0 0 0 0 0 0 0 1 0 0 0
45 9788 0 0 0 0 0 0 0 0 1 0 0
46 10312 0 0 0 0 0 0 0 0 0 1 0
47 10105 0 0 0 0 0 0 0 0 0 0 1
48 9863 0 0 0 0 0 0 0 0 0 0 0
49 9656 1 0 0 0 0 0 0 0 0 0 0
50 9295 0 1 0 0 0 0 0 0 0 0 0
51 9946 0 0 1 0 0 0 0 0 0 0 0
52 9701 0 0 0 1 0 0 0 0 0 0 0
53 9049 0 0 0 0 1 0 0 0 0 0 0
54 10190 0 0 0 0 0 1 0 0 0 0 0
55 9706 0 0 0 0 0 0 1 0 0 0 0
56 9765 0 0 0 0 0 0 0 1 0 0 0
57 9893 0 0 0 0 0 0 0 0 1 0 0
58 9994 0 0 0 0 0 0 0 0 0 1 0
59 10433 0 0 0 0 0 0 0 0 0 0 1
60 10073 0 0 0 0 0 0 0 0 0 0 0
61 10112 1 0 0 0 0 0 0 0 0 0 0
62 9266 0 1 0 0 0 0 0 0 0 0 0
63 9820 0 0 1 0 0 0 0 0 0 0 0
64 10097 0 0 0 1 0 0 0 0 0 0 0
65 9115 0 0 0 0 1 0 0 0 0 0 0
66 10411 0 0 0 0 0 1 0 0 0 0 0
67 9678 0 0 0 0 0 0 1 0 0 0 0
68 10408 0 0 0 0 0 0 0 1 0 0 0
69 10153 0 0 0 0 0 0 0 0 1 0 0
70 10368 0 0 0 0 0 0 0 0 0 1 0
71 10581 0 0 0 0 0 0 0 0 0 0 1
72 10597 0 0 0 0 0 0 0 0 0 0 0
73 10680 1 0 0 0 0 0 0 0 0 0 0
74 9738 0 1 0 0 0 0 0 0 0 0 0
75 9556 0 0 1 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
9793.000 52.571 -679.571 -320.857 -79.333 -956.833
M6 M7 M8 M9 M10 M11
9.667 -364.667 -185.333 -230.000 345.167 223.500
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-734.00 -244.87 -32.43 239.75 834.43
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9793.000 158.732 61.695 < 2e-16 ***
M1 52.571 216.316 0.243 0.80877
M2 -679.571 216.316 -3.142 0.00256 **
M3 -320.857 216.316 -1.483 0.14299
M4 -79.333 224.482 -0.353 0.72496
M5 -956.833 224.482 -4.262 6.89e-05 ***
M6 9.667 224.482 0.043 0.96579
M7 -364.667 224.482 -1.624 0.10927
M8 -185.333 224.482 -0.826 0.41214
M9 -230.000 224.482 -1.025 0.30948
M10 345.167 224.482 1.538 0.12915
M11 223.500 224.482 0.996 0.32324
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 388.8 on 63 degrees of freedom
Multiple R-squared: 0.4914, Adjusted R-squared: 0.4026
F-statistic: 5.535 on 11 and 63 DF, p-value: 4.088e-06
> 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.0005841378 0.0011682756 0.9994159
[2,] 0.0048917514 0.0097835028 0.9951082
[3,] 0.0013759063 0.0027518125 0.9986241
[4,] 0.0003671796 0.0007343592 0.9996328
[5,] 0.0131169821 0.0262339641 0.9868830
[6,] 0.0519488335 0.1038976670 0.9480512
[7,] 0.0937355692 0.1874711385 0.9062644
[8,] 0.0565177568 0.1130355137 0.9434822
[9,] 0.0430969787 0.0861939575 0.9569030
[10,] 0.0306359947 0.0612719894 0.9693640
[11,] 0.0186369494 0.0372738988 0.9813631
[12,] 0.0238930934 0.0477861868 0.9761069
[13,] 0.0160348376 0.0320696752 0.9839652
[14,] 0.0090664455 0.0181328910 0.9909336
[15,] 0.0059140508 0.0118281015 0.9940859
[16,] 0.0151609100 0.0303218200 0.9848391
[17,] 0.0100072375 0.0200144751 0.9899928
[18,] 0.0187076574 0.0374153148 0.9812923
[19,] 0.0174508826 0.0349017653 0.9825491
[20,] 0.0112068402 0.0224136803 0.9887932
[21,] 0.0235050174 0.0470100348 0.9764950
[22,] 0.0282712760 0.0565425521 0.9717287
[23,] 0.0395815470 0.0791630940 0.9604185
[24,] 0.0521789414 0.1043578827 0.9478211
[25,] 0.0549462828 0.1098925655 0.9450537
[26,] 0.0387798598 0.0775597196 0.9612201
[27,] 0.0372331067 0.0744662133 0.9627669
[28,] 0.0438800226 0.0877600452 0.9561200
[29,] 0.0417868515 0.0835737030 0.9582131
[30,] 0.1868823197 0.3737646394 0.8131177
[31,] 0.2391066297 0.4782132593 0.7608934
[32,] 0.1990238204 0.3980476408 0.8009762
[33,] 0.2324067979 0.4648135958 0.7675932
[34,] 0.2929953964 0.5859907929 0.7070046
[35,] 0.4917224980 0.9834449960 0.5082775
[36,] 0.4745148402 0.9490296804 0.5254852
[37,] 0.5378655667 0.9242688666 0.4621344
[38,] 0.5188789329 0.9622421342 0.4811211
[39,] 0.4435903866 0.8871807732 0.5564096
[40,] 0.4390328967 0.8780657935 0.5609671
[41,] 0.3615475789 0.7230951578 0.6384524
[42,] 0.4897162814 0.9794325628 0.5102837
[43,] 0.4483292879 0.8966585759 0.5516707
[44,] 0.4047978580 0.8095957161 0.5952021
[45,] 0.3256415952 0.6512831903 0.6743584
[46,] 0.3725795560 0.7451591119 0.6274204
> postscript(file="/var/fisher/rcomp/tmp/1dp7r1356129886.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/2wjsz1356129886.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/37gx71356129886.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/4697o1356129886.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/51tdb1356129886.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 7
-145.57143 -32.42857 -388.14286 29.33333 -249.16667 -71.66667 134.66667
8 9 10 11 12 13 14
390.33333 -126.00000 -100.16667 -98.50000 -541.00000 -108.57143 -78.42857
15 16 17 18 19 20 21
-339.14286 -226.66667 -136.16667 -175.66667 -481.33333 -324.66667 -734.00000
22 23 24 25 26 27 28
-191.16667 -388.50000 -475.00000 -240.57143 -473.42857 -258.14286 -146.66667
29 30 31 32 33 34 35
-289.16667 -617.66667 41.66667 -484.66667 -285.00000 31.83333 -582.50000
36 37 38 39 40 41 42
-138.00000 -416.57143 -374.42857 79.85714 -26.66667 182.83333 -130.66667
43 44 45 46 47 48 49
-222.33333 -538.66667 225.00000 173.83333 88.50000 70.00000 -189.57143
50 51 52 53 54 55 56
181.57143 473.85714 -12.66667 212.83333 387.33333 277.66667 157.33333
57 58 59 60 61 62 63
330.00000 -144.16667 416.50000 280.00000 266.42857 152.57143 347.85714
64 65 66 67 68 69 70
383.33333 278.83333 608.33333 249.66667 800.33333 590.00000 229.83333
71 72 73 74 75
564.50000 804.00000 834.42857 624.57143 83.85714
> postscript(file="/var/fisher/rcomp/tmp/6umrh1356129886.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 -145.57143 NA
1 -32.42857 -145.57143
2 -388.14286 -32.42857
3 29.33333 -388.14286
4 -249.16667 29.33333
5 -71.66667 -249.16667
6 134.66667 -71.66667
7 390.33333 134.66667
8 -126.00000 390.33333
9 -100.16667 -126.00000
10 -98.50000 -100.16667
11 -541.00000 -98.50000
12 -108.57143 -541.00000
13 -78.42857 -108.57143
14 -339.14286 -78.42857
15 -226.66667 -339.14286
16 -136.16667 -226.66667
17 -175.66667 -136.16667
18 -481.33333 -175.66667
19 -324.66667 -481.33333
20 -734.00000 -324.66667
21 -191.16667 -734.00000
22 -388.50000 -191.16667
23 -475.00000 -388.50000
24 -240.57143 -475.00000
25 -473.42857 -240.57143
26 -258.14286 -473.42857
27 -146.66667 -258.14286
28 -289.16667 -146.66667
29 -617.66667 -289.16667
30 41.66667 -617.66667
31 -484.66667 41.66667
32 -285.00000 -484.66667
33 31.83333 -285.00000
34 -582.50000 31.83333
35 -138.00000 -582.50000
36 -416.57143 -138.00000
37 -374.42857 -416.57143
38 79.85714 -374.42857
39 -26.66667 79.85714
40 182.83333 -26.66667
41 -130.66667 182.83333
42 -222.33333 -130.66667
43 -538.66667 -222.33333
44 225.00000 -538.66667
45 173.83333 225.00000
46 88.50000 173.83333
47 70.00000 88.50000
48 -189.57143 70.00000
49 181.57143 -189.57143
50 473.85714 181.57143
51 -12.66667 473.85714
52 212.83333 -12.66667
53 387.33333 212.83333
54 277.66667 387.33333
55 157.33333 277.66667
56 330.00000 157.33333
57 -144.16667 330.00000
58 416.50000 -144.16667
59 280.00000 416.50000
60 266.42857 280.00000
61 152.57143 266.42857
62 347.85714 152.57143
63 383.33333 347.85714
64 278.83333 383.33333
65 608.33333 278.83333
66 249.66667 608.33333
67 800.33333 249.66667
68 590.00000 800.33333
69 229.83333 590.00000
70 564.50000 229.83333
71 804.00000 564.50000
72 834.42857 804.00000
73 624.57143 834.42857
74 83.85714 624.57143
75 NA 83.85714
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -32.42857 -145.57143
[2,] -388.14286 -32.42857
[3,] 29.33333 -388.14286
[4,] -249.16667 29.33333
[5,] -71.66667 -249.16667
[6,] 134.66667 -71.66667
[7,] 390.33333 134.66667
[8,] -126.00000 390.33333
[9,] -100.16667 -126.00000
[10,] -98.50000 -100.16667
[11,] -541.00000 -98.50000
[12,] -108.57143 -541.00000
[13,] -78.42857 -108.57143
[14,] -339.14286 -78.42857
[15,] -226.66667 -339.14286
[16,] -136.16667 -226.66667
[17,] -175.66667 -136.16667
[18,] -481.33333 -175.66667
[19,] -324.66667 -481.33333
[20,] -734.00000 -324.66667
[21,] -191.16667 -734.00000
[22,] -388.50000 -191.16667
[23,] -475.00000 -388.50000
[24,] -240.57143 -475.00000
[25,] -473.42857 -240.57143
[26,] -258.14286 -473.42857
[27,] -146.66667 -258.14286
[28,] -289.16667 -146.66667
[29,] -617.66667 -289.16667
[30,] 41.66667 -617.66667
[31,] -484.66667 41.66667
[32,] -285.00000 -484.66667
[33,] 31.83333 -285.00000
[34,] -582.50000 31.83333
[35,] -138.00000 -582.50000
[36,] -416.57143 -138.00000
[37,] -374.42857 -416.57143
[38,] 79.85714 -374.42857
[39,] -26.66667 79.85714
[40,] 182.83333 -26.66667
[41,] -130.66667 182.83333
[42,] -222.33333 -130.66667
[43,] -538.66667 -222.33333
[44,] 225.00000 -538.66667
[45,] 173.83333 225.00000
[46,] 88.50000 173.83333
[47,] 70.00000 88.50000
[48,] -189.57143 70.00000
[49,] 181.57143 -189.57143
[50,] 473.85714 181.57143
[51,] -12.66667 473.85714
[52,] 212.83333 -12.66667
[53,] 387.33333 212.83333
[54,] 277.66667 387.33333
[55,] 157.33333 277.66667
[56,] 330.00000 157.33333
[57,] -144.16667 330.00000
[58,] 416.50000 -144.16667
[59,] 280.00000 416.50000
[60,] 266.42857 280.00000
[61,] 152.57143 266.42857
[62,] 347.85714 152.57143
[63,] 383.33333 347.85714
[64,] 278.83333 383.33333
[65,] 608.33333 278.83333
[66,] 249.66667 608.33333
[67,] 800.33333 249.66667
[68,] 590.00000 800.33333
[69,] 229.83333 590.00000
[70,] 564.50000 229.83333
[71,] 804.00000 564.50000
[72,] 834.42857 804.00000
[73,] 624.57143 834.42857
[74,] 83.85714 624.57143
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -32.42857 -145.57143
2 -388.14286 -32.42857
3 29.33333 -388.14286
4 -249.16667 29.33333
5 -71.66667 -249.16667
6 134.66667 -71.66667
7 390.33333 134.66667
8 -126.00000 390.33333
9 -100.16667 -126.00000
10 -98.50000 -100.16667
11 -541.00000 -98.50000
12 -108.57143 -541.00000
13 -78.42857 -108.57143
14 -339.14286 -78.42857
15 -226.66667 -339.14286
16 -136.16667 -226.66667
17 -175.66667 -136.16667
18 -481.33333 -175.66667
19 -324.66667 -481.33333
20 -734.00000 -324.66667
21 -191.16667 -734.00000
22 -388.50000 -191.16667
23 -475.00000 -388.50000
24 -240.57143 -475.00000
25 -473.42857 -240.57143
26 -258.14286 -473.42857
27 -146.66667 -258.14286
28 -289.16667 -146.66667
29 -617.66667 -289.16667
30 41.66667 -617.66667
31 -484.66667 41.66667
32 -285.00000 -484.66667
33 31.83333 -285.00000
34 -582.50000 31.83333
35 -138.00000 -582.50000
36 -416.57143 -138.00000
37 -374.42857 -416.57143
38 79.85714 -374.42857
39 -26.66667 79.85714
40 182.83333 -26.66667
41 -130.66667 182.83333
42 -222.33333 -130.66667
43 -538.66667 -222.33333
44 225.00000 -538.66667
45 173.83333 225.00000
46 88.50000 173.83333
47 70.00000 88.50000
48 -189.57143 70.00000
49 181.57143 -189.57143
50 473.85714 181.57143
51 -12.66667 473.85714
52 212.83333 -12.66667
53 387.33333 212.83333
54 277.66667 387.33333
55 157.33333 277.66667
56 330.00000 157.33333
57 -144.16667 330.00000
58 416.50000 -144.16667
59 280.00000 416.50000
60 266.42857 280.00000
61 152.57143 266.42857
62 347.85714 152.57143
63 383.33333 347.85714
64 278.83333 383.33333
65 608.33333 278.83333
66 249.66667 608.33333
67 800.33333 249.66667
68 590.00000 800.33333
69 229.83333 590.00000
70 564.50000 229.83333
71 804.00000 564.50000
72 834.42857 804.00000
73 624.57143 834.42857
74 83.85714 624.57143
> 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/782sl1356129886.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/82opk1356129886.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/96ect1356129886.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/10dok21356129886.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/1192qw1356129886.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/1231251356129886.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/13ytel1356129886.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/14wfjt1356129886.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/15rl1t1356129887.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/16okt41356129887.tab")
+ }
>
> try(system("convert tmp/1dp7r1356129886.ps tmp/1dp7r1356129886.png",intern=TRUE))
character(0)
> try(system("convert tmp/2wjsz1356129886.ps tmp/2wjsz1356129886.png",intern=TRUE))
character(0)
> try(system("convert tmp/37gx71356129886.ps tmp/37gx71356129886.png",intern=TRUE))
character(0)
> try(system("convert tmp/4697o1356129886.ps tmp/4697o1356129886.png",intern=TRUE))
character(0)
> try(system("convert tmp/51tdb1356129886.ps tmp/51tdb1356129886.png",intern=TRUE))
character(0)
> try(system("convert tmp/6umrh1356129886.ps tmp/6umrh1356129886.png",intern=TRUE))
character(0)
> try(system("convert tmp/782sl1356129886.ps tmp/782sl1356129886.png",intern=TRUE))
character(0)
> try(system("convert tmp/82opk1356129886.ps tmp/82opk1356129886.png",intern=TRUE))
character(0)
> try(system("convert tmp/96ect1356129886.ps tmp/96ect1356129886.png",intern=TRUE))
character(0)
> try(system("convert tmp/10dok21356129886.ps tmp/10dok21356129886.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
5.985 1.688 7.671