R version 2.9.0 (2009-04-17)
Copyright (C) 2009 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.
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 = '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
Garnalen M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 3010 1 0 0 0 0 0 0 0 0 0 0
2 2910 0 1 0 0 0 0 0 0 0 0 0
3 3840 0 0 1 0 0 0 0 0 0 0 0
4 3580 0 0 0 1 0 0 0 0 0 0 0
5 3140 0 0 0 0 1 0 0 0 0 0 0
6 3550 0 0 0 0 0 1 0 0 0 0 0
7 3250 0 0 0 0 0 0 1 0 0 0 0
8 2820 0 0 0 0 0 0 0 1 0 0 0
9 2260 0 0 0 0 0 0 0 0 1 0 0
10 2060 0 0 0 0 0 0 0 0 0 1 0
11 2120 0 0 0 0 0 0 0 0 0 0 1
12 2210 0 0 0 0 0 0 0 0 0 0 0
13 2190 1 0 0 0 0 0 0 0 0 0 0
14 2180 0 1 0 0 0 0 0 0 0 0 0
15 2350 0 0 1 0 0 0 0 0 0 0 0
16 2440 0 0 0 1 0 0 0 0 0 0 0
17 2370 0 0 0 0 1 0 0 0 0 0 0
18 2440 0 0 0 0 0 1 0 0 0 0 0
19 2610 0 0 0 0 0 0 1 0 0 0 0
20 3040 0 0 0 0 0 0 0 1 0 0 0
21 3190 0 0 0 0 0 0 0 0 1 0 0
22 3120 0 0 0 0 0 0 0 0 0 1 0
23 3170 0 0 0 0 0 0 0 0 0 0 1
24 3600 0 0 0 0 0 0 0 0 0 0 0
25 3420 1 0 0 0 0 0 0 0 0 0 0
26 3650 0 1 0 0 0 0 0 0 0 0 0
27 4180 0 0 1 0 0 0 0 0 0 0 0
28 2960 0 0 0 1 0 0 0 0 0 0 0
29 2710 0 0 0 0 1 0 0 0 0 0 0
30 2950 0 0 0 0 0 1 0 0 0 0 0
31 3030 0 0 0 0 0 0 1 0 0 0 0
32 3770 0 0 0 0 0 0 0 1 0 0 0
33 4740 0 0 0 0 0 0 0 0 1 0 0
34 4450 0 0 0 0 0 0 0 0 0 1 0
35 5550 0 0 0 0 0 0 0 0 0 0 1
36 5580 0 0 0 0 0 0 0 0 0 0 0
37 5890 1 0 0 0 0 0 0 0 0 0 0
38 7480 0 1 0 0 0 0 0 0 0 0 0
39 10450 0 0 1 0 0 0 0 0 0 0 0
40 6360 0 0 0 1 0 0 0 0 0 0 0
41 6710 0 0 0 0 1 0 0 0 0 0 0
42 6200 0 0 0 0 0 1 0 0 0 0 0
43 4490 0 0 0 0 0 0 1 0 0 0 0
44 3480 0 0 0 0 0 0 0 1 0 0 0
45 2520 0 0 0 0 0 0 0 0 1 0 0
46 1920 0 0 0 0 0 0 0 0 0 1 0
47 2010 0 0 0 0 0 0 0 0 0 0 1
48 1950 0 0 0 0 0 0 0 0 0 0 0
49 2240 1 0 0 0 0 0 0 0 0 0 0
50 2370 0 1 0 0 0 0 0 0 0 0 0
51 2840 0 0 1 0 0 0 0 0 0 0 0
52 2700 0 0 0 1 0 0 0 0 0 0 0
53 2980 0 0 0 0 1 0 0 0 0 0 0
54 3290 0 0 0 0 0 1 0 0 0 0 0
55 3300 0 0 0 0 0 0 1 0 0 0 0
56 3000 0 0 0 0 0 0 0 1 0 0 0
57 2330 0 0 0 0 0 0 0 0 1 0 0
58 2190 0 0 0 0 0 0 0 0 0 1 0
59 1970 0 0 0 0 0 0 0 0 0 0 1
60 2170 0 0 0 0 0 0 0 0 0 0 0
61 2830 1 0 0 0 0 0 0 0 0 0 0
62 3190 0 1 0 0 0 0 0 0 0 0 0
63 3550 0 0 1 0 0 0 0 0 0 0 0
64 3240 0 0 0 1 0 0 0 0 0 0 0
65 3450 0 0 0 0 1 0 0 0 0 0 0
66 3570 0 0 0 0 0 1 0 0 0 0 0
67 3230 0 0 0 0 0 0 1 0 0 0 0
68 3260 0 0 0 0 0 0 0 1 0 0 0
69 2700 0 0 0 0 0 0 0 0 1 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
3102.0 161.3 528.0 1433.0 444.7 458.0
M6 M7 M8 M9 M10 M11
564.7 216.3 126.3 -145.3 -354.0 -138.0
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2185.0 -844.0 -376.7 156.7 5915.0
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3102.0 686.9 4.516 3.22e-05 ***
M1 161.3 930.1 0.173 0.863
M2 528.0 930.1 0.568 0.572
M3 1433.0 930.1 1.541 0.129
M4 444.7 930.1 0.478 0.634
M5 458.0 930.1 0.492 0.624
M6 564.7 930.1 0.607 0.546
M7 216.3 930.1 0.233 0.817
M8 126.3 930.1 0.136 0.892
M9 -145.3 930.1 -0.156 0.876
M10 -354.0 971.4 -0.364 0.717
M11 -138.0 971.4 -0.142 0.888
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1536 on 57 degrees of freedom
Multiple R-squared: 0.09297, Adjusted R-squared: -0.08207
F-statistic: 0.5311 on 11 and 57 DF, p-value: 0.874
> 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,] 1.061702e-01 2.123403e-01 8.938298e-01
[2,] 6.721709e-02 1.344342e-01 9.327829e-01
[3,] 3.357158e-02 6.714315e-02 9.664284e-01
[4,] 2.182096e-02 4.364193e-02 9.781790e-01
[5,] 9.950736e-03 1.990147e-02 9.900493e-01
[6,] 3.669109e-03 7.338218e-03 9.963309e-01
[7,] 1.952735e-03 3.905469e-03 9.980473e-01
[8,] 1.169762e-03 2.339524e-03 9.988302e-01
[9,] 6.831371e-04 1.366274e-03 9.993169e-01
[10,] 5.677253e-04 1.135451e-03 9.994323e-01
[11,] 2.952008e-04 5.904016e-04 9.997048e-01
[12,] 2.082011e-04 4.164022e-04 9.997918e-01
[13,] 1.497833e-04 2.995665e-04 9.998502e-01
[14,] 5.599967e-05 1.119993e-04 9.999440e-01
[15,] 2.172789e-05 4.345578e-05 9.999783e-01
[16,] 7.876892e-06 1.575378e-05 9.999921e-01
[17,] 2.547709e-06 5.095417e-06 9.999975e-01
[18,] 1.261467e-06 2.522935e-06 9.999987e-01
[19,] 5.641432e-06 1.128286e-05 9.999944e-01
[20,] 1.330160e-05 2.660321e-05 9.999867e-01
[21,] 1.909662e-04 3.819323e-04 9.998090e-01
[22,] 9.834911e-04 1.966982e-03 9.990165e-01
[23,] 5.535587e-03 1.107117e-02 9.944644e-01
[24,] 9.894544e-02 1.978909e-01 9.010546e-01
[25,] 9.704041e-01 5.919180e-02 2.959590e-02
[26,] 9.959526e-01 8.094813e-03 4.047407e-03
[27,] 9.999547e-01 9.067000e-05 4.533500e-05
[28,] 1.000000e+00 6.160151e-08 3.080075e-08
[29,] 1.000000e+00 6.212186e-09 3.106093e-09
[30,] 1.000000e+00 2.653009e-08 1.326505e-08
[31,] 9.999999e-01 1.624497e-07 8.122484e-08
[32,] 9.999996e-01 7.826094e-07 3.913047e-07
[33,] 9.999979e-01 4.162757e-06 2.081379e-06
[34,] 9.999905e-01 1.905399e-05 9.526993e-06
[35,] 9.999744e-01 5.121898e-05 2.560949e-05
[36,] 9.999684e-01 6.312706e-05 3.156353e-05
[37,] 9.999581e-01 8.381920e-05 4.190960e-05
[38,] 9.999007e-01 1.985744e-04 9.928721e-05
[39,] 9.997307e-01 5.386956e-04 2.693478e-04
[40,] 9.983203e-01 3.359367e-03 1.679683e-03
> postscript(file="/var/www/html/rcomp/tmp/1x34i1292945330.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/rcomp/tmp/2x34i1292945330.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/rcomp/tmp/3x34i1292945330.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/rcomp/tmp/4qcl31292945330.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/rcomp/tmp/5qcl31292945330.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
-253.33333 -720.00000 -695.00000 33.33333 -420.00000 -116.66667
7 8 9 10 11 12
-68.33333 -408.33333 -696.66667 -688.00000 -844.00000 -892.00000
13 14 15 16 17 18
-1073.33333 -1450.00000 -2185.00000 -1106.66667 -1190.00000 -1226.66667
19 20 21 22 23 24
-708.33333 -188.33333 233.33333 372.00000 206.00000 498.00000
25 26 27 28 29 30
156.66667 20.00000 -355.00000 -586.66667 -850.00000 -716.66667
31 32 33 34 35 36
-288.33333 541.66667 1783.33333 1702.00000 2586.00000 2478.00000
37 38 39 40 41 42
2626.66667 3850.00000 5915.00000 2813.33333 3150.00000 2533.33333
43 44 45 46 47 48
1171.66667 251.66667 -436.66667 -828.00000 -954.00000 -1152.00000
49 50 51 52 53 54
-1023.33333 -1260.00000 -1695.00000 -846.66667 -580.00000 -376.66667
55 56 57 58 59 60
-18.33333 -228.33333 -626.66667 -558.00000 -994.00000 -932.00000
61 62 63 64 65 66
-433.33333 -440.00000 -985.00000 -306.66667 -110.00000 -96.66667
67 68 69
-88.33333 31.66667 -256.66667
> postscript(file="/var/www/html/rcomp/tmp/6qcl31292945330.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 -253.33333 NA
1 -720.00000 -253.33333
2 -695.00000 -720.00000
3 33.33333 -695.00000
4 -420.00000 33.33333
5 -116.66667 -420.00000
6 -68.33333 -116.66667
7 -408.33333 -68.33333
8 -696.66667 -408.33333
9 -688.00000 -696.66667
10 -844.00000 -688.00000
11 -892.00000 -844.00000
12 -1073.33333 -892.00000
13 -1450.00000 -1073.33333
14 -2185.00000 -1450.00000
15 -1106.66667 -2185.00000
16 -1190.00000 -1106.66667
17 -1226.66667 -1190.00000
18 -708.33333 -1226.66667
19 -188.33333 -708.33333
20 233.33333 -188.33333
21 372.00000 233.33333
22 206.00000 372.00000
23 498.00000 206.00000
24 156.66667 498.00000
25 20.00000 156.66667
26 -355.00000 20.00000
27 -586.66667 -355.00000
28 -850.00000 -586.66667
29 -716.66667 -850.00000
30 -288.33333 -716.66667
31 541.66667 -288.33333
32 1783.33333 541.66667
33 1702.00000 1783.33333
34 2586.00000 1702.00000
35 2478.00000 2586.00000
36 2626.66667 2478.00000
37 3850.00000 2626.66667
38 5915.00000 3850.00000
39 2813.33333 5915.00000
40 3150.00000 2813.33333
41 2533.33333 3150.00000
42 1171.66667 2533.33333
43 251.66667 1171.66667
44 -436.66667 251.66667
45 -828.00000 -436.66667
46 -954.00000 -828.00000
47 -1152.00000 -954.00000
48 -1023.33333 -1152.00000
49 -1260.00000 -1023.33333
50 -1695.00000 -1260.00000
51 -846.66667 -1695.00000
52 -580.00000 -846.66667
53 -376.66667 -580.00000
54 -18.33333 -376.66667
55 -228.33333 -18.33333
56 -626.66667 -228.33333
57 -558.00000 -626.66667
58 -994.00000 -558.00000
59 -932.00000 -994.00000
60 -433.33333 -932.00000
61 -440.00000 -433.33333
62 -985.00000 -440.00000
63 -306.66667 -985.00000
64 -110.00000 -306.66667
65 -96.66667 -110.00000
66 -88.33333 -96.66667
67 31.66667 -88.33333
68 -256.66667 31.66667
69 NA -256.66667
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -720.00000 -253.33333
[2,] -695.00000 -720.00000
[3,] 33.33333 -695.00000
[4,] -420.00000 33.33333
[5,] -116.66667 -420.00000
[6,] -68.33333 -116.66667
[7,] -408.33333 -68.33333
[8,] -696.66667 -408.33333
[9,] -688.00000 -696.66667
[10,] -844.00000 -688.00000
[11,] -892.00000 -844.00000
[12,] -1073.33333 -892.00000
[13,] -1450.00000 -1073.33333
[14,] -2185.00000 -1450.00000
[15,] -1106.66667 -2185.00000
[16,] -1190.00000 -1106.66667
[17,] -1226.66667 -1190.00000
[18,] -708.33333 -1226.66667
[19,] -188.33333 -708.33333
[20,] 233.33333 -188.33333
[21,] 372.00000 233.33333
[22,] 206.00000 372.00000
[23,] 498.00000 206.00000
[24,] 156.66667 498.00000
[25,] 20.00000 156.66667
[26,] -355.00000 20.00000
[27,] -586.66667 -355.00000
[28,] -850.00000 -586.66667
[29,] -716.66667 -850.00000
[30,] -288.33333 -716.66667
[31,] 541.66667 -288.33333
[32,] 1783.33333 541.66667
[33,] 1702.00000 1783.33333
[34,] 2586.00000 1702.00000
[35,] 2478.00000 2586.00000
[36,] 2626.66667 2478.00000
[37,] 3850.00000 2626.66667
[38,] 5915.00000 3850.00000
[39,] 2813.33333 5915.00000
[40,] 3150.00000 2813.33333
[41,] 2533.33333 3150.00000
[42,] 1171.66667 2533.33333
[43,] 251.66667 1171.66667
[44,] -436.66667 251.66667
[45,] -828.00000 -436.66667
[46,] -954.00000 -828.00000
[47,] -1152.00000 -954.00000
[48,] -1023.33333 -1152.00000
[49,] -1260.00000 -1023.33333
[50,] -1695.00000 -1260.00000
[51,] -846.66667 -1695.00000
[52,] -580.00000 -846.66667
[53,] -376.66667 -580.00000
[54,] -18.33333 -376.66667
[55,] -228.33333 -18.33333
[56,] -626.66667 -228.33333
[57,] -558.00000 -626.66667
[58,] -994.00000 -558.00000
[59,] -932.00000 -994.00000
[60,] -433.33333 -932.00000
[61,] -440.00000 -433.33333
[62,] -985.00000 -440.00000
[63,] -306.66667 -985.00000
[64,] -110.00000 -306.66667
[65,] -96.66667 -110.00000
[66,] -88.33333 -96.66667
[67,] 31.66667 -88.33333
[68,] -256.66667 31.66667
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -720.00000 -253.33333
2 -695.00000 -720.00000
3 33.33333 -695.00000
4 -420.00000 33.33333
5 -116.66667 -420.00000
6 -68.33333 -116.66667
7 -408.33333 -68.33333
8 -696.66667 -408.33333
9 -688.00000 -696.66667
10 -844.00000 -688.00000
11 -892.00000 -844.00000
12 -1073.33333 -892.00000
13 -1450.00000 -1073.33333
14 -2185.00000 -1450.00000
15 -1106.66667 -2185.00000
16 -1190.00000 -1106.66667
17 -1226.66667 -1190.00000
18 -708.33333 -1226.66667
19 -188.33333 -708.33333
20 233.33333 -188.33333
21 372.00000 233.33333
22 206.00000 372.00000
23 498.00000 206.00000
24 156.66667 498.00000
25 20.00000 156.66667
26 -355.00000 20.00000
27 -586.66667 -355.00000
28 -850.00000 -586.66667
29 -716.66667 -850.00000
30 -288.33333 -716.66667
31 541.66667 -288.33333
32 1783.33333 541.66667
33 1702.00000 1783.33333
34 2586.00000 1702.00000
35 2478.00000 2586.00000
36 2626.66667 2478.00000
37 3850.00000 2626.66667
38 5915.00000 3850.00000
39 2813.33333 5915.00000
40 3150.00000 2813.33333
41 2533.33333 3150.00000
42 1171.66667 2533.33333
43 251.66667 1171.66667
44 -436.66667 251.66667
45 -828.00000 -436.66667
46 -954.00000 -828.00000
47 -1152.00000 -954.00000
48 -1023.33333 -1152.00000
49 -1260.00000 -1023.33333
50 -1695.00000 -1260.00000
51 -846.66667 -1695.00000
52 -580.00000 -846.66667
53 -376.66667 -580.00000
54 -18.33333 -376.66667
55 -228.33333 -18.33333
56 -626.66667 -228.33333
57 -558.00000 -626.66667
58 -994.00000 -558.00000
59 -932.00000 -994.00000
60 -433.33333 -932.00000
61 -440.00000 -433.33333
62 -985.00000 -440.00000
63 -306.66667 -985.00000
64 -110.00000 -306.66667
65 -96.66667 -110.00000
66 -88.33333 -96.66667
67 31.66667 -88.33333
68 -256.66667 31.66667
> 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/rcomp/tmp/7i32o1292945330.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/rcomp/tmp/8bdk91292945330.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/rcomp/tmp/9bdk91292945330.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/rcomp/tmp/10mmju1292945330.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/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/rcomp/tmp/11uqoo1292945330.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/rcomp/tmp/12iwh21292945330.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/rcomp/tmp/137xww1292945330.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/rcomp/tmp/14axc21292945330.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/rcomp/tmp/153pc51292945330.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/rcomp/tmp/16hz9e1292945330.tab")
+ }
>
> try(system("convert tmp/1x34i1292945330.ps tmp/1x34i1292945330.png",intern=TRUE))
character(0)
> try(system("convert tmp/2x34i1292945330.ps tmp/2x34i1292945330.png",intern=TRUE))
character(0)
> try(system("convert tmp/3x34i1292945330.ps tmp/3x34i1292945330.png",intern=TRUE))
character(0)
> try(system("convert tmp/4qcl31292945330.ps tmp/4qcl31292945330.png",intern=TRUE))
character(0)
> try(system("convert tmp/5qcl31292945330.ps tmp/5qcl31292945330.png",intern=TRUE))
character(0)
> try(system("convert tmp/6qcl31292945330.ps tmp/6qcl31292945330.png",intern=TRUE))
character(0)
> try(system("convert tmp/7i32o1292945330.ps tmp/7i32o1292945330.png",intern=TRUE))
character(0)
> try(system("convert tmp/8bdk91292945330.ps tmp/8bdk91292945330.png",intern=TRUE))
character(0)
> try(system("convert tmp/9bdk91292945330.ps tmp/9bdk91292945330.png",intern=TRUE))
character(0)
> try(system("convert tmp/10mmju1292945330.ps tmp/10mmju1292945330.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.614 1.627 9.638