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(3010,2910,3840,3580,3140,3550,3250,2820,2260,2060,2120,2210,2190,2180,2350,2440,2370,2440,2610,3040,3190,3120,3170,3600,3420,3650,4180,2960,2710,2950,3030,3770,4740,4450,5550,5580,5890,7480,10450,6360,6710,6200,4490,3480,2520,1920,2010,1950,2240,2370,2840,2700,2980,3290,3300,3000,2330,2190,1970,2170,2830,3190,3550,3240,3450,3570,3230,3260,2700),dim=c(1,69),dimnames=list(c('Garnalen'),1:69))
> y <- array(NA,dim=c(1,69),dimnames=list(c('Garnalen'),1:69))
> 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 = '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
Garnalen M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 3010 1 0 0 0 0 0 0 0 0 0 0 1
2 2910 0 1 0 0 0 0 0 0 0 0 0 2
3 3840 0 0 1 0 0 0 0 0 0 0 0 3
4 3580 0 0 0 1 0 0 0 0 0 0 0 4
5 3140 0 0 0 0 1 0 0 0 0 0 0 5
6 3550 0 0 0 0 0 1 0 0 0 0 0 6
7 3250 0 0 0 0 0 0 1 0 0 0 0 7
8 2820 0 0 0 0 0 0 0 1 0 0 0 8
9 2260 0 0 0 0 0 0 0 0 1 0 0 9
10 2060 0 0 0 0 0 0 0 0 0 1 0 10
11 2120 0 0 0 0 0 0 0 0 0 0 1 11
12 2210 0 0 0 0 0 0 0 0 0 0 0 12
13 2190 1 0 0 0 0 0 0 0 0 0 0 13
14 2180 0 1 0 0 0 0 0 0 0 0 0 14
15 2350 0 0 1 0 0 0 0 0 0 0 0 15
16 2440 0 0 0 1 0 0 0 0 0 0 0 16
17 2370 0 0 0 0 1 0 0 0 0 0 0 17
18 2440 0 0 0 0 0 1 0 0 0 0 0 18
19 2610 0 0 0 0 0 0 1 0 0 0 0 19
20 3040 0 0 0 0 0 0 0 1 0 0 0 20
21 3190 0 0 0 0 0 0 0 0 1 0 0 21
22 3120 0 0 0 0 0 0 0 0 0 1 0 22
23 3170 0 0 0 0 0 0 0 0 0 0 1 23
24 3600 0 0 0 0 0 0 0 0 0 0 0 24
25 3420 1 0 0 0 0 0 0 0 0 0 0 25
26 3650 0 1 0 0 0 0 0 0 0 0 0 26
27 4180 0 0 1 0 0 0 0 0 0 0 0 27
28 2960 0 0 0 1 0 0 0 0 0 0 0 28
29 2710 0 0 0 0 1 0 0 0 0 0 0 29
30 2950 0 0 0 0 0 1 0 0 0 0 0 30
31 3030 0 0 0 0 0 0 1 0 0 0 0 31
32 3770 0 0 0 0 0 0 0 1 0 0 0 32
33 4740 0 0 0 0 0 0 0 0 1 0 0 33
34 4450 0 0 0 0 0 0 0 0 0 1 0 34
35 5550 0 0 0 0 0 0 0 0 0 0 1 35
36 5580 0 0 0 0 0 0 0 0 0 0 0 36
37 5890 1 0 0 0 0 0 0 0 0 0 0 37
38 7480 0 1 0 0 0 0 0 0 0 0 0 38
39 10450 0 0 1 0 0 0 0 0 0 0 0 39
40 6360 0 0 0 1 0 0 0 0 0 0 0 40
41 6710 0 0 0 0 1 0 0 0 0 0 0 41
42 6200 0 0 0 0 0 1 0 0 0 0 0 42
43 4490 0 0 0 0 0 0 1 0 0 0 0 43
44 3480 0 0 0 0 0 0 0 1 0 0 0 44
45 2520 0 0 0 0 0 0 0 0 1 0 0 45
46 1920 0 0 0 0 0 0 0 0 0 1 0 46
47 2010 0 0 0 0 0 0 0 0 0 0 1 47
48 1950 0 0 0 0 0 0 0 0 0 0 0 48
49 2240 1 0 0 0 0 0 0 0 0 0 0 49
50 2370 0 1 0 0 0 0 0 0 0 0 0 50
51 2840 0 0 1 0 0 0 0 0 0 0 0 51
52 2700 0 0 0 1 0 0 0 0 0 0 0 52
53 2980 0 0 0 0 1 0 0 0 0 0 0 53
54 3290 0 0 0 0 0 1 0 0 0 0 0 54
55 3300 0 0 0 0 0 0 1 0 0 0 0 55
56 3000 0 0 0 0 0 0 0 1 0 0 0 56
57 2330 0 0 0 0 0 0 0 0 1 0 0 57
58 2190 0 0 0 0 0 0 0 0 0 1 0 58
59 1970 0 0 0 0 0 0 0 0 0 0 1 59
60 2170 0 0 0 0 0 0 0 0 0 0 0 60
61 2830 1 0 0 0 0 0 0 0 0 0 0 61
62 3190 0 1 0 0 0 0 0 0 0 0 0 62
63 3550 0 0 1 0 0 0 0 0 0 0 0 63
64 3240 0 0 0 1 0 0 0 0 0 0 0 64
65 3450 0 0 0 0 1 0 0 0 0 0 0 65
66 3570 0 0 0 0 0 1 0 0 0 0 0 66
67 3230 0 0 0 0 0 0 1 0 0 0 0 67
68 3260 0 0 0 0 0 0 0 1 0 0 0 68
69 2700 0 0 0 0 0 0 0 0 1 0 0 69
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
2910.560 187.922 549.271 1448.953 455.302 463.318
M6 M7 M8 M9 M10 M11
564.667 211.016 115.698 -161.287 -343.364 -132.682
t
5.318
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2089.3 -764.4 -416.2 192.9 5883.1
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2910.560 769.511 3.782 0.00038 ***
M1 187.922 936.855 0.201 0.84175
M2 549.271 936.430 0.587 0.55986
M3 1448.953 936.099 1.548 0.12729
M4 455.302 935.863 0.487 0.62851
M5 463.318 935.721 0.495 0.62244
M6 564.667 935.674 0.603 0.54862
M7 211.016 935.721 0.226 0.82240
M8 115.698 935.863 0.124 0.90205
M9 -161.287 936.099 -0.172 0.86383
M10 -343.364 977.461 -0.351 0.72669
M11 -132.682 977.325 -0.136 0.89250
t 5.318 9.404 0.565 0.57400
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1545 on 56 degrees of freedom
Multiple R-squared: 0.09812, Adjusted R-squared: -0.09514
F-statistic: 0.5077 on 12 and 56 DF, p-value: 0.901
> 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,] 4.258205e-03 8.516411e-03 9.957418e-01
[2,] 7.230016e-04 1.446003e-03 9.992770e-01
[3,] 9.812039e-05 1.962408e-04 9.999019e-01
[4,] 2.328538e-05 4.657077e-05 9.999767e-01
[5,] 1.494581e-04 2.989162e-04 9.998505e-01
[6,] 9.467384e-04 1.893477e-03 9.990533e-01
[7,] 1.636303e-03 3.272607e-03 9.983637e-01
[8,] 1.673476e-03 3.346951e-03 9.983265e-01
[9,] 1.957651e-03 3.915303e-03 9.980423e-01
[10,] 1.294047e-03 2.588094e-03 9.987060e-01
[11,] 9.985320e-04 1.997064e-03 9.990015e-01
[12,] 8.005075e-04 1.601015e-03 9.991995e-01
[13,] 4.796992e-04 9.593984e-04 9.995203e-01
[14,] 3.976694e-04 7.953388e-04 9.996023e-01
[15,] 3.831009e-04 7.662019e-04 9.996169e-01
[16,] 3.307151e-04 6.614303e-04 9.996693e-01
[17,] 2.584018e-04 5.168036e-04 9.997416e-01
[18,] 5.083223e-04 1.016645e-03 9.994917e-01
[19,] 5.536559e-04 1.107312e-03 9.994463e-01
[20,] 2.142192e-03 4.284384e-03 9.978578e-01
[21,] 3.767690e-03 7.535380e-03 9.962323e-01
[22,] 5.885952e-03 1.177190e-02 9.941140e-01
[23,] 4.222430e-02 8.444859e-02 9.577757e-01
[24,] 8.852558e-01 2.294884e-01 1.147442e-01
[25,] 9.486542e-01 1.026916e-01 5.134579e-02
[26,] 9.961165e-01 7.767096e-03 3.883548e-03
[27,] 9.999869e-01 2.628442e-05 1.314221e-05
[28,] 9.999998e-01 3.934422e-07 1.967211e-07
[29,] 1.000000e+00 9.685332e-08 4.842666e-08
[30,] 1.000000e+00 5.393759e-08 2.696880e-08
[31,] 9.999999e-01 2.121112e-07 1.060556e-07
[32,] 9.999999e-01 2.915227e-07 1.457613e-07
[33,] 9.999994e-01 1.218944e-06 6.094719e-07
[34,] 9.999966e-01 6.747160e-06 3.373580e-06
[35,] 9.999930e-01 1.390415e-05 6.952076e-06
[36,] 9.999830e-01 3.400819e-05 1.700409e-05
[37,] 9.999003e-01 1.993768e-04 9.968840e-05
[38,] 9.993084e-01 1.383187e-03 6.915935e-04
> postscript(file="/var/www/html/freestat/rcomp/tmp/17w2m1292946564.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/27w2m1292946564.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/37w2m1292946564.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/4injp1292946564.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/5injp1292946564.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 = 69
Frequency = 1
1 2 3 4 5 6
-93.80000 -560.46667 -535.46667 192.86667 -260.46667 42.86667
7 8 9 10 11 12
91.20000 -248.80000 -537.13333 -560.37333 -716.37333 -764.37333
13 14 15 16 17 18
-977.61333 -1354.28000 -2089.28000 -1010.94667 -1094.28000 -1130.94667
19 20 21 22 23 24
-612.61333 -92.61333 329.05333 435.81333 269.81333 561.81333
25 26 27 28 29 30
188.57333 51.90667 -323.09333 -554.76000 -818.09333 -684.76000
31 32 33 34 35 36
-256.42667 573.57333 1815.24000 1702.00000 2586.00000 2478.00000
37 38 39 40 41 42
2594.76000 3818.09333 5883.09333 2781.42667 3118.09333 2501.42667
43 44 45 46 47 48
1139.76000 219.76000 -468.57333 -891.81333 -1017.81333 -1215.81333
49 50 51 52 53 54
-1119.05333 -1355.72000 -1790.72000 -942.38667 -675.72000 -472.38667
55 56 57 58 59 60
-114.05333 -324.05333 -722.38667 -685.62667 -1121.62667 -1059.62667
61 62 63 64 65 66
-592.86667 -599.53333 -1144.53333 -466.20000 -269.53333 -256.20000
67 68 69
-247.86667 -127.86667 -416.20000
> postscript(file="/var/www/html/freestat/rcomp/tmp/6injp1292946564.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 = 69
Frequency = 1
lag(myerror, k = 1) myerror
0 -93.80000 NA
1 -560.46667 -93.80000
2 -535.46667 -560.46667
3 192.86667 -535.46667
4 -260.46667 192.86667
5 42.86667 -260.46667
6 91.20000 42.86667
7 -248.80000 91.20000
8 -537.13333 -248.80000
9 -560.37333 -537.13333
10 -716.37333 -560.37333
11 -764.37333 -716.37333
12 -977.61333 -764.37333
13 -1354.28000 -977.61333
14 -2089.28000 -1354.28000
15 -1010.94667 -2089.28000
16 -1094.28000 -1010.94667
17 -1130.94667 -1094.28000
18 -612.61333 -1130.94667
19 -92.61333 -612.61333
20 329.05333 -92.61333
21 435.81333 329.05333
22 269.81333 435.81333
23 561.81333 269.81333
24 188.57333 561.81333
25 51.90667 188.57333
26 -323.09333 51.90667
27 -554.76000 -323.09333
28 -818.09333 -554.76000
29 -684.76000 -818.09333
30 -256.42667 -684.76000
31 573.57333 -256.42667
32 1815.24000 573.57333
33 1702.00000 1815.24000
34 2586.00000 1702.00000
35 2478.00000 2586.00000
36 2594.76000 2478.00000
37 3818.09333 2594.76000
38 5883.09333 3818.09333
39 2781.42667 5883.09333
40 3118.09333 2781.42667
41 2501.42667 3118.09333
42 1139.76000 2501.42667
43 219.76000 1139.76000
44 -468.57333 219.76000
45 -891.81333 -468.57333
46 -1017.81333 -891.81333
47 -1215.81333 -1017.81333
48 -1119.05333 -1215.81333
49 -1355.72000 -1119.05333
50 -1790.72000 -1355.72000
51 -942.38667 -1790.72000
52 -675.72000 -942.38667
53 -472.38667 -675.72000
54 -114.05333 -472.38667
55 -324.05333 -114.05333
56 -722.38667 -324.05333
57 -685.62667 -722.38667
58 -1121.62667 -685.62667
59 -1059.62667 -1121.62667
60 -592.86667 -1059.62667
61 -599.53333 -592.86667
62 -1144.53333 -599.53333
63 -466.20000 -1144.53333
64 -269.53333 -466.20000
65 -256.20000 -269.53333
66 -247.86667 -256.20000
67 -127.86667 -247.86667
68 -416.20000 -127.86667
69 NA -416.20000
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -560.46667 -93.80000
[2,] -535.46667 -560.46667
[3,] 192.86667 -535.46667
[4,] -260.46667 192.86667
[5,] 42.86667 -260.46667
[6,] 91.20000 42.86667
[7,] -248.80000 91.20000
[8,] -537.13333 -248.80000
[9,] -560.37333 -537.13333
[10,] -716.37333 -560.37333
[11,] -764.37333 -716.37333
[12,] -977.61333 -764.37333
[13,] -1354.28000 -977.61333
[14,] -2089.28000 -1354.28000
[15,] -1010.94667 -2089.28000
[16,] -1094.28000 -1010.94667
[17,] -1130.94667 -1094.28000
[18,] -612.61333 -1130.94667
[19,] -92.61333 -612.61333
[20,] 329.05333 -92.61333
[21,] 435.81333 329.05333
[22,] 269.81333 435.81333
[23,] 561.81333 269.81333
[24,] 188.57333 561.81333
[25,] 51.90667 188.57333
[26,] -323.09333 51.90667
[27,] -554.76000 -323.09333
[28,] -818.09333 -554.76000
[29,] -684.76000 -818.09333
[30,] -256.42667 -684.76000
[31,] 573.57333 -256.42667
[32,] 1815.24000 573.57333
[33,] 1702.00000 1815.24000
[34,] 2586.00000 1702.00000
[35,] 2478.00000 2586.00000
[36,] 2594.76000 2478.00000
[37,] 3818.09333 2594.76000
[38,] 5883.09333 3818.09333
[39,] 2781.42667 5883.09333
[40,] 3118.09333 2781.42667
[41,] 2501.42667 3118.09333
[42,] 1139.76000 2501.42667
[43,] 219.76000 1139.76000
[44,] -468.57333 219.76000
[45,] -891.81333 -468.57333
[46,] -1017.81333 -891.81333
[47,] -1215.81333 -1017.81333
[48,] -1119.05333 -1215.81333
[49,] -1355.72000 -1119.05333
[50,] -1790.72000 -1355.72000
[51,] -942.38667 -1790.72000
[52,] -675.72000 -942.38667
[53,] -472.38667 -675.72000
[54,] -114.05333 -472.38667
[55,] -324.05333 -114.05333
[56,] -722.38667 -324.05333
[57,] -685.62667 -722.38667
[58,] -1121.62667 -685.62667
[59,] -1059.62667 -1121.62667
[60,] -592.86667 -1059.62667
[61,] -599.53333 -592.86667
[62,] -1144.53333 -599.53333
[63,] -466.20000 -1144.53333
[64,] -269.53333 -466.20000
[65,] -256.20000 -269.53333
[66,] -247.86667 -256.20000
[67,] -127.86667 -247.86667
[68,] -416.20000 -127.86667
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -560.46667 -93.80000
2 -535.46667 -560.46667
3 192.86667 -535.46667
4 -260.46667 192.86667
5 42.86667 -260.46667
6 91.20000 42.86667
7 -248.80000 91.20000
8 -537.13333 -248.80000
9 -560.37333 -537.13333
10 -716.37333 -560.37333
11 -764.37333 -716.37333
12 -977.61333 -764.37333
13 -1354.28000 -977.61333
14 -2089.28000 -1354.28000
15 -1010.94667 -2089.28000
16 -1094.28000 -1010.94667
17 -1130.94667 -1094.28000
18 -612.61333 -1130.94667
19 -92.61333 -612.61333
20 329.05333 -92.61333
21 435.81333 329.05333
22 269.81333 435.81333
23 561.81333 269.81333
24 188.57333 561.81333
25 51.90667 188.57333
26 -323.09333 51.90667
27 -554.76000 -323.09333
28 -818.09333 -554.76000
29 -684.76000 -818.09333
30 -256.42667 -684.76000
31 573.57333 -256.42667
32 1815.24000 573.57333
33 1702.00000 1815.24000
34 2586.00000 1702.00000
35 2478.00000 2586.00000
36 2594.76000 2478.00000
37 3818.09333 2594.76000
38 5883.09333 3818.09333
39 2781.42667 5883.09333
40 3118.09333 2781.42667
41 2501.42667 3118.09333
42 1139.76000 2501.42667
43 219.76000 1139.76000
44 -468.57333 219.76000
45 -891.81333 -468.57333
46 -1017.81333 -891.81333
47 -1215.81333 -1017.81333
48 -1119.05333 -1215.81333
49 -1355.72000 -1119.05333
50 -1790.72000 -1355.72000
51 -942.38667 -1790.72000
52 -675.72000 -942.38667
53 -472.38667 -675.72000
54 -114.05333 -472.38667
55 -324.05333 -114.05333
56 -722.38667 -324.05333
57 -685.62667 -722.38667
58 -1121.62667 -685.62667
59 -1059.62667 -1121.62667
60 -592.86667 -1059.62667
61 -599.53333 -592.86667
62 -1144.53333 -599.53333
63 -466.20000 -1144.53333
64 -269.53333 -466.20000
65 -256.20000 -269.53333
66 -247.86667 -256.20000
67 -127.86667 -247.86667
68 -416.20000 -127.86667
> 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/7awi91292946564.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/835id1292946564.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/935id1292946564.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/1035id1292946564.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/11hgjv1292946565.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/12kyzj1292946565.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/139hwv1292946565.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/1429eg1292946565.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/155rc41292946565.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/161jau1292946565.tab")
+ }
>
> try(system("convert tmp/17w2m1292946564.ps tmp/17w2m1292946564.png",intern=TRUE))
character(0)
> try(system("convert tmp/27w2m1292946564.ps tmp/27w2m1292946564.png",intern=TRUE))
character(0)
> try(system("convert tmp/37w2m1292946564.ps tmp/37w2m1292946564.png",intern=TRUE))
character(0)
> try(system("convert tmp/4injp1292946564.ps tmp/4injp1292946564.png",intern=TRUE))
character(0)
> try(system("convert tmp/5injp1292946564.ps tmp/5injp1292946564.png",intern=TRUE))
character(0)
> try(system("convert tmp/6injp1292946564.ps tmp/6injp1292946564.png",intern=TRUE))
character(0)
> try(system("convert tmp/7awi91292946564.ps tmp/7awi91292946564.png",intern=TRUE))
character(0)
> try(system("convert tmp/835id1292946564.ps tmp/835id1292946564.png",intern=TRUE))
character(0)
> try(system("convert tmp/935id1292946564.ps tmp/935id1292946564.png",intern=TRUE))
character(0)
> try(system("convert tmp/1035id1292946564.ps tmp/1035id1292946564.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.904 2.473 4.339