R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(9700,0,9081,0,9084,0,9743,0,8587,0,9731,0,9563,0,9998,0,9437,0,10038,0,9918,0,9252,0,9737,0,9035,0,9133,0,9487,0,8700,0,9627,0,8947,0,9283,0,8829,0,9947,0,9628,0,9318,0,9605,0,8640,0,9214,0,9567,0,8547,0,9185,0,9470,0,9123,0,9278,0,10170,0,9434,0,9655,0,9429,0,8739,0,9552,0,9687,1,9019,1,9672,1,9206,1,9069,1,9788,1,10312,1,10105,1,9863,1,9656,1,9295,1,9946,1,9701,1,9049,1,10190,1,9706,1,9765,1,9893,1,9994,1,10433,1,10073,1,10112,1,9266,1,9820,1,10097,1,9115,1,10411,1,9678,1,10408,1,10153,1,10368,1,10581,1,10597,1,10680,1,9738,1,9556,1),dim=c(2,75),dimnames=list(c('births','difference'),1:75))
> y <- array(NA,dim=c(2,75),dimnames=list(c('births','difference'),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'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
births difference M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 9700 0 1 0 0 0 0 0 0 0 0 0 0
2 9081 0 0 1 0 0 0 0 0 0 0 0 0
3 9084 0 0 0 1 0 0 0 0 0 0 0 0
4 9743 0 0 0 0 1 0 0 0 0 0 0 0
5 8587 0 0 0 0 0 1 0 0 0 0 0 0
6 9731 0 0 0 0 0 0 1 0 0 0 0 0
7 9563 0 0 0 0 0 0 0 1 0 0 0 0
8 9998 0 0 0 0 0 0 0 0 1 0 0 0
9 9437 0 0 0 0 0 0 0 0 0 1 0 0
10 10038 0 0 0 0 0 0 0 0 0 0 1 0
11 9918 0 0 0 0 0 0 0 0 0 0 0 1
12 9252 0 0 0 0 0 0 0 0 0 0 0 0
13 9737 0 1 0 0 0 0 0 0 0 0 0 0
14 9035 0 0 1 0 0 0 0 0 0 0 0 0
15 9133 0 0 0 1 0 0 0 0 0 0 0 0
16 9487 0 0 0 0 1 0 0 0 0 0 0 0
17 8700 0 0 0 0 0 1 0 0 0 0 0 0
18 9627 0 0 0 0 0 0 1 0 0 0 0 0
19 8947 0 0 0 0 0 0 0 1 0 0 0 0
20 9283 0 0 0 0 0 0 0 0 1 0 0 0
21 8829 0 0 0 0 0 0 0 0 0 1 0 0
22 9947 0 0 0 0 0 0 0 0 0 0 1 0
23 9628 0 0 0 0 0 0 0 0 0 0 0 1
24 9318 0 0 0 0 0 0 0 0 0 0 0 0
25 9605 0 1 0 0 0 0 0 0 0 0 0 0
26 8640 0 0 1 0 0 0 0 0 0 0 0 0
27 9214 0 0 0 1 0 0 0 0 0 0 0 0
28 9567 0 0 0 0 1 0 0 0 0 0 0 0
29 8547 0 0 0 0 0 1 0 0 0 0 0 0
30 9185 0 0 0 0 0 0 1 0 0 0 0 0
31 9470 0 0 0 0 0 0 0 1 0 0 0 0
32 9123 0 0 0 0 0 0 0 0 1 0 0 0
33 9278 0 0 0 0 0 0 0 0 0 1 0 0
34 10170 0 0 0 0 0 0 0 0 0 0 1 0
35 9434 0 0 0 0 0 0 0 0 0 0 0 1
36 9655 0 0 0 0 0 0 0 0 0 0 0 0
37 9429 0 1 0 0 0 0 0 0 0 0 0 0
38 8739 0 0 1 0 0 0 0 0 0 0 0 0
39 9552 0 0 0 1 0 0 0 0 0 0 0 0
40 9687 1 0 0 0 1 0 0 0 0 0 0 0
41 9019 1 0 0 0 0 1 0 0 0 0 0 0
42 9672 1 0 0 0 0 0 1 0 0 0 0 0
43 9206 1 0 0 0 0 0 0 1 0 0 0 0
44 9069 1 0 0 0 0 0 0 0 1 0 0 0
45 9788 1 0 0 0 0 0 0 0 0 1 0 0
46 10312 1 0 0 0 0 0 0 0 0 0 1 0
47 10105 1 0 0 0 0 0 0 0 0 0 0 1
48 9863 1 0 0 0 0 0 0 0 0 0 0 0
49 9656 1 1 0 0 0 0 0 0 0 0 0 0
50 9295 1 0 1 0 0 0 0 0 0 0 0 0
51 9946 1 0 0 1 0 0 0 0 0 0 0 0
52 9701 1 0 0 0 1 0 0 0 0 0 0 0
53 9049 1 0 0 0 0 1 0 0 0 0 0 0
54 10190 1 0 0 0 0 0 1 0 0 0 0 0
55 9706 1 0 0 0 0 0 0 1 0 0 0 0
56 9765 1 0 0 0 0 0 0 0 1 0 0 0
57 9893 1 0 0 0 0 0 0 0 0 1 0 0
58 9994 1 0 0 0 0 0 0 0 0 0 1 0
59 10433 1 0 0 0 0 0 0 0 0 0 0 1
60 10073 1 0 0 0 0 0 0 0 0 0 0 0
61 10112 1 1 0 0 0 0 0 0 0 0 0 0
62 9266 1 0 1 0 0 0 0 0 0 0 0 0
63 9820 1 0 0 1 0 0 0 0 0 0 0 0
64 10097 1 0 0 0 1 0 0 0 0 0 0 0
65 9115 1 0 0 0 0 1 0 0 0 0 0 0
66 10411 1 0 0 0 0 0 1 0 0 0 0 0
67 9678 1 0 0 0 0 0 0 1 0 0 0 0
68 10408 1 0 0 0 0 0 0 0 1 0 0 0
69 10153 1 0 0 0 0 0 0 0 0 1 0 0
70 10368 1 0 0 0 0 0 0 0 0 0 1 0
71 10581 1 0 0 0 0 0 0 0 0 0 0 1
72 10597 1 0 0 0 0 0 0 0 0 0 0 0
73 10680 1 1 0 0 0 0 0 0 0 0 0 0
74 9738 1 0 1 0 0 0 0 0 0 0 0 0
75 9556 1 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) difference M1 M2 M3 M4
9551.324 483.352 87.097 -645.046 -286.332 -79.333
M5 M6 M7 M8 M9 M10
-956.833 9.667 -364.667 -185.333 -230.000 345.167
M11
223.500
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-780.34 -169.48 -7.49 142.42 632.01
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9551.324 122.523 77.955 < 2e-16 ***
difference 483.352 66.870 7.228 8.65e-10 ***
M1 87.097 160.704 0.542 0.589783
M2 -645.046 160.704 -4.014 0.000164 ***
M3 -286.332 160.704 -1.782 0.079691 .
M4 -79.333 166.697 -0.476 0.635809
M5 -956.833 166.697 -5.740 3.05e-07 ***
M6 9.667 166.697 0.058 0.953944
M7 -364.667 166.697 -2.188 0.032478 *
M8 -185.333 166.697 -1.112 0.270518
M9 -230.000 166.697 -1.380 0.172620
M10 345.167 166.697 2.071 0.042564 *
M11 223.500 166.697 1.341 0.184892
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 288.7 on 62 degrees of freedom
Multiple R-squared: 0.724, Adjusted R-squared: 0.6706
F-statistic: 13.55 on 12 and 62 DF, p-value: 3.705e-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.062886409 0.12577282 0.9371136
[2,] 0.025253838 0.05050768 0.9747462
[3,] 0.009709587 0.01941917 0.9902904
[4,] 0.148579320 0.29715864 0.8514207
[5,] 0.375555161 0.75111032 0.6244448
[6,] 0.508910415 0.98217917 0.4910896
[7,] 0.406094554 0.81218911 0.5939054
[8,] 0.346973812 0.69394762 0.6530262
[9,] 0.272384862 0.54476972 0.7276151
[10,] 0.201453921 0.40290784 0.7985461
[11,] 0.220629085 0.44125817 0.7793709
[12,] 0.162010428 0.32402086 0.8379896
[13,] 0.117628050 0.23525610 0.8823719
[14,] 0.080264637 0.16052927 0.9197354
[15,] 0.116285925 0.23257185 0.8837141
[16,] 0.107301041 0.21460208 0.8926990
[17,] 0.134188896 0.26837779 0.8658111
[18,] 0.099058934 0.19811787 0.9009411
[19,] 0.094790070 0.18958014 0.9052099
[20,] 0.095133733 0.19026747 0.9048663
[21,] 0.086249306 0.17249861 0.9137507
[22,] 0.071270264 0.14254053 0.9287297
[23,] 0.058964925 0.11792985 0.9410351
[24,] 0.057185784 0.11437157 0.9428142
[25,] 0.039314598 0.07862920 0.9606854
[26,] 0.028282817 0.05656563 0.9717172
[27,] 0.031417791 0.06283558 0.9685822
[28,] 0.034791411 0.06958282 0.9652086
[29,] 0.176210808 0.35242162 0.8237892
[30,] 0.201813127 0.40362625 0.7981869
[31,] 0.154250354 0.30850071 0.8457496
[32,] 0.158510708 0.31702142 0.8414893
[33,] 0.179424777 0.35884955 0.8205752
[34,] 0.350709780 0.70141956 0.6492902
[35,] 0.303553555 0.60710711 0.6964464
[36,] 0.310627321 0.62125464 0.6893727
[37,] 0.295266014 0.59053203 0.7047340
[38,] 0.219145808 0.43829162 0.7808542
[39,] 0.194213920 0.38842784 0.8057861
[40,] 0.132094428 0.26418886 0.8679056
[41,] 0.207467474 0.41493495 0.7925325
[42,] 0.165151569 0.33030314 0.8348484
[43,] 0.151190529 0.30238106 0.8488095
[44,] 0.095786838 0.19157368 0.9042132
> postscript(file="/var/www/html/freestat/rcomp/tmp/1383z1291915546.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/2dhk21291915546.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/3dhk21291915546.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/4dhk21291915546.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/5dhk21291915546.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
61.579639 174.722496 -180.991790 271.009579 -7.490421 170.009579
7 8 9 10 11 12
376.342912 632.009579 115.676245 141.509579 143.176245 -299.323755
13 14 15 16 17 18
98.579639 128.722496 -131.991790 15.009579 105.509579 66.009579
19 20 21 22 23 24
-239.657088 -82.990421 -492.323755 50.509579 -146.823755 -233.323755
25 26 27 28 29 30
-33.420361 -266.277504 -50.991790 95.009579 -47.490421 -375.990421
31 32 33 34 35 36
283.342912 -242.990421 -43.323755 273.509579 -340.823755 103.676245
37 38 39 40 41 42
-209.420361 -167.277504 287.008210 -268.342912 -58.842912 -372.342912
43 44 45 46 47 48
-464.009579 -780.342912 -16.676245 -67.842912 -153.176245 -171.676245
49 50 51 52 53 54
-465.772852 -94.629995 197.655720 -254.342912 -28.842912 145.657088
55 56 57 58 59 60
35.990421 -84.342912 88.323755 -385.842912 174.823755 38.323755
61 62 63 64 65 66
-9.772852 -123.629995 71.655720 141.657088 37.157088 366.657088
67 68 69 70 71 72
7.990421 558.657088 348.323755 -11.842912 322.823755 562.323755
73 74 75
558.227148 348.370005 -192.344280
> postscript(file="/var/www/html/freestat/rcomp/tmp/6682n1291915546.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 61.579639 NA
1 174.722496 61.579639
2 -180.991790 174.722496
3 271.009579 -180.991790
4 -7.490421 271.009579
5 170.009579 -7.490421
6 376.342912 170.009579
7 632.009579 376.342912
8 115.676245 632.009579
9 141.509579 115.676245
10 143.176245 141.509579
11 -299.323755 143.176245
12 98.579639 -299.323755
13 128.722496 98.579639
14 -131.991790 128.722496
15 15.009579 -131.991790
16 105.509579 15.009579
17 66.009579 105.509579
18 -239.657088 66.009579
19 -82.990421 -239.657088
20 -492.323755 -82.990421
21 50.509579 -492.323755
22 -146.823755 50.509579
23 -233.323755 -146.823755
24 -33.420361 -233.323755
25 -266.277504 -33.420361
26 -50.991790 -266.277504
27 95.009579 -50.991790
28 -47.490421 95.009579
29 -375.990421 -47.490421
30 283.342912 -375.990421
31 -242.990421 283.342912
32 -43.323755 -242.990421
33 273.509579 -43.323755
34 -340.823755 273.509579
35 103.676245 -340.823755
36 -209.420361 103.676245
37 -167.277504 -209.420361
38 287.008210 -167.277504
39 -268.342912 287.008210
40 -58.842912 -268.342912
41 -372.342912 -58.842912
42 -464.009579 -372.342912
43 -780.342912 -464.009579
44 -16.676245 -780.342912
45 -67.842912 -16.676245
46 -153.176245 -67.842912
47 -171.676245 -153.176245
48 -465.772852 -171.676245
49 -94.629995 -465.772852
50 197.655720 -94.629995
51 -254.342912 197.655720
52 -28.842912 -254.342912
53 145.657088 -28.842912
54 35.990421 145.657088
55 -84.342912 35.990421
56 88.323755 -84.342912
57 -385.842912 88.323755
58 174.823755 -385.842912
59 38.323755 174.823755
60 -9.772852 38.323755
61 -123.629995 -9.772852
62 71.655720 -123.629995
63 141.657088 71.655720
64 37.157088 141.657088
65 366.657088 37.157088
66 7.990421 366.657088
67 558.657088 7.990421
68 348.323755 558.657088
69 -11.842912 348.323755
70 322.823755 -11.842912
71 562.323755 322.823755
72 558.227148 562.323755
73 348.370005 558.227148
74 -192.344280 348.370005
75 NA -192.344280
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 174.722496 61.579639
[2,] -180.991790 174.722496
[3,] 271.009579 -180.991790
[4,] -7.490421 271.009579
[5,] 170.009579 -7.490421
[6,] 376.342912 170.009579
[7,] 632.009579 376.342912
[8,] 115.676245 632.009579
[9,] 141.509579 115.676245
[10,] 143.176245 141.509579
[11,] -299.323755 143.176245
[12,] 98.579639 -299.323755
[13,] 128.722496 98.579639
[14,] -131.991790 128.722496
[15,] 15.009579 -131.991790
[16,] 105.509579 15.009579
[17,] 66.009579 105.509579
[18,] -239.657088 66.009579
[19,] -82.990421 -239.657088
[20,] -492.323755 -82.990421
[21,] 50.509579 -492.323755
[22,] -146.823755 50.509579
[23,] -233.323755 -146.823755
[24,] -33.420361 -233.323755
[25,] -266.277504 -33.420361
[26,] -50.991790 -266.277504
[27,] 95.009579 -50.991790
[28,] -47.490421 95.009579
[29,] -375.990421 -47.490421
[30,] 283.342912 -375.990421
[31,] -242.990421 283.342912
[32,] -43.323755 -242.990421
[33,] 273.509579 -43.323755
[34,] -340.823755 273.509579
[35,] 103.676245 -340.823755
[36,] -209.420361 103.676245
[37,] -167.277504 -209.420361
[38,] 287.008210 -167.277504
[39,] -268.342912 287.008210
[40,] -58.842912 -268.342912
[41,] -372.342912 -58.842912
[42,] -464.009579 -372.342912
[43,] -780.342912 -464.009579
[44,] -16.676245 -780.342912
[45,] -67.842912 -16.676245
[46,] -153.176245 -67.842912
[47,] -171.676245 -153.176245
[48,] -465.772852 -171.676245
[49,] -94.629995 -465.772852
[50,] 197.655720 -94.629995
[51,] -254.342912 197.655720
[52,] -28.842912 -254.342912
[53,] 145.657088 -28.842912
[54,] 35.990421 145.657088
[55,] -84.342912 35.990421
[56,] 88.323755 -84.342912
[57,] -385.842912 88.323755
[58,] 174.823755 -385.842912
[59,] 38.323755 174.823755
[60,] -9.772852 38.323755
[61,] -123.629995 -9.772852
[62,] 71.655720 -123.629995
[63,] 141.657088 71.655720
[64,] 37.157088 141.657088
[65,] 366.657088 37.157088
[66,] 7.990421 366.657088
[67,] 558.657088 7.990421
[68,] 348.323755 558.657088
[69,] -11.842912 348.323755
[70,] 322.823755 -11.842912
[71,] 562.323755 322.823755
[72,] 558.227148 562.323755
[73,] 348.370005 558.227148
[74,] -192.344280 348.370005
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 174.722496 61.579639
2 -180.991790 174.722496
3 271.009579 -180.991790
4 -7.490421 271.009579
5 170.009579 -7.490421
6 376.342912 170.009579
7 632.009579 376.342912
8 115.676245 632.009579
9 141.509579 115.676245
10 143.176245 141.509579
11 -299.323755 143.176245
12 98.579639 -299.323755
13 128.722496 98.579639
14 -131.991790 128.722496
15 15.009579 -131.991790
16 105.509579 15.009579
17 66.009579 105.509579
18 -239.657088 66.009579
19 -82.990421 -239.657088
20 -492.323755 -82.990421
21 50.509579 -492.323755
22 -146.823755 50.509579
23 -233.323755 -146.823755
24 -33.420361 -233.323755
25 -266.277504 -33.420361
26 -50.991790 -266.277504
27 95.009579 -50.991790
28 -47.490421 95.009579
29 -375.990421 -47.490421
30 283.342912 -375.990421
31 -242.990421 283.342912
32 -43.323755 -242.990421
33 273.509579 -43.323755
34 -340.823755 273.509579
35 103.676245 -340.823755
36 -209.420361 103.676245
37 -167.277504 -209.420361
38 287.008210 -167.277504
39 -268.342912 287.008210
40 -58.842912 -268.342912
41 -372.342912 -58.842912
42 -464.009579 -372.342912
43 -780.342912 -464.009579
44 -16.676245 -780.342912
45 -67.842912 -16.676245
46 -153.176245 -67.842912
47 -171.676245 -153.176245
48 -465.772852 -171.676245
49 -94.629995 -465.772852
50 197.655720 -94.629995
51 -254.342912 197.655720
52 -28.842912 -254.342912
53 145.657088 -28.842912
54 35.990421 145.657088
55 -84.342912 35.990421
56 88.323755 -84.342912
57 -385.842912 88.323755
58 174.823755 -385.842912
59 38.323755 174.823755
60 -9.772852 38.323755
61 -123.629995 -9.772852
62 71.655720 -123.629995
63 141.657088 71.655720
64 37.157088 141.657088
65 366.657088 37.157088
66 7.990421 366.657088
67 558.657088 7.990421
68 348.323755 558.657088
69 -11.842912 348.323755
70 322.823755 -11.842912
71 562.323755 322.823755
72 558.227148 562.323755
73 348.370005 558.227148
74 -192.344280 348.370005
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/7zi1q1291915546.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/8zi1q1291915546.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/9rr0t1291915546.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/freestat/rcomp/tmp/10rr0t1291915546.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/11vrzh1291915546.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/12gaxm1291915546.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/13sejn1291915546.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/14g2b11291915546.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/1513a71291915546.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/freestat/rcomp/tmp/164lqv1291915546.tab")
+ }
>
> try(system("convert tmp/1383z1291915546.ps tmp/1383z1291915546.png",intern=TRUE))
character(0)
> try(system("convert tmp/2dhk21291915546.ps tmp/2dhk21291915546.png",intern=TRUE))
character(0)
> try(system("convert tmp/3dhk21291915546.ps tmp/3dhk21291915546.png",intern=TRUE))
character(0)
> try(system("convert tmp/4dhk21291915546.ps tmp/4dhk21291915546.png",intern=TRUE))
character(0)
> try(system("convert tmp/5dhk21291915546.ps tmp/5dhk21291915546.png",intern=TRUE))
character(0)
> try(system("convert tmp/6682n1291915546.ps tmp/6682n1291915546.png",intern=TRUE))
character(0)
> try(system("convert tmp/7zi1q1291915546.ps tmp/7zi1q1291915546.png",intern=TRUE))
character(0)
> try(system("convert tmp/8zi1q1291915546.ps tmp/8zi1q1291915546.png",intern=TRUE))
character(0)
> try(system("convert tmp/9rr0t1291915546.ps tmp/9rr0t1291915546.png",intern=TRUE))
character(0)
> try(system("convert tmp/10rr0t1291915546.ps tmp/10rr0t1291915546.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.029 2.519 4.354