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
+ ,2590
+ ,11290
+ ,4700
+ ,44.51
+ ,2910
+ ,2080
+ ,11620
+ ,4960
+ ,45.48
+ ,3840
+ ,2640
+ ,10790
+ ,4880
+ ,53.1
+ ,3580
+ ,3000
+ ,8380
+ ,4090
+ ,51.88
+ ,3140
+ ,2350
+ ,9370
+ ,3450
+ ,48.65
+ ,3550
+ ,2220
+ ,10090
+ ,3020
+ ,54.35
+ ,3250
+ ,2540
+ ,11130
+ ,3070
+ ,57.52
+ ,2820
+ ,2700
+ ,10530
+ ,3720
+ ,63.98
+ ,2260
+ ,2580
+ ,10490
+ ,3750
+ ,62.91
+ ,2060
+ ,2420
+ ,10570
+ ,3910
+ ,58.54
+ ,2120
+ ,2090
+ ,11170
+ ,4120
+ ,55.24
+ ,2210
+ ,2000
+ ,11610
+ ,4780
+ ,56.86
+ ,2190
+ ,1860
+ ,10920
+ ,3070
+ ,62.99
+ ,2180
+ ,1980
+ ,11570
+ ,4100
+ ,60.21
+ ,2350
+ ,2690
+ ,12960
+ ,3900
+ ,62.06
+ ,2440
+ ,3040
+ ,11190
+ ,3020
+ ,70.26
+ ,2370
+ ,2450
+ ,11920
+ ,3220
+ ,69.78
+ ,2440
+ ,2650
+ ,14930
+ ,4030
+ ,68.56
+ ,2610
+ ,2710
+ ,14520
+ ,4210
+ ,73.67
+ ,3040
+ ,3230
+ ,12970
+ ,4510
+ ,73.23
+ ,3190
+ ,3160
+ ,13870
+ ,4320
+ ,61.96
+ ,3120
+ ,3040
+ ,13250
+ ,3890
+ ,57.81
+ ,3170
+ ,2630
+ ,12760
+ ,7280
+ ,58.76
+ ,3600
+ ,2730
+ ,14050
+ ,9640
+ ,62.47
+ ,3420
+ ,2830
+ ,14660
+ ,5680
+ ,53.68
+ ,3650
+ ,2320
+ ,15010
+ ,6320
+ ,57.56
+ ,4180
+ ,2410
+ ,15020
+ ,5820
+ ,62.05
+ ,2960
+ ,3080
+ ,13090
+ ,4890
+ ,67.49
+ ,2710
+ ,2260
+ ,13190
+ ,3320
+ ,67.21
+ ,2950
+ ,2300
+ ,11390
+ ,2930
+ ,71.05
+ ,3030
+ ,3600
+ ,10110
+ ,3530
+ ,76.93
+ ,3770
+ ,3380
+ ,8240
+ ,3690
+ ,70.76
+ ,4740
+ ,3670
+ ,7920
+ ,3750
+ ,77.17
+ ,4450
+ ,3040
+ ,7700
+ ,3330
+ ,82.34
+ ,5550
+ ,2840
+ ,7920
+ ,4790
+ ,92.41
+ ,5580
+ ,2810
+ ,8130
+ ,5990
+ ,90.93
+ ,5890
+ ,2980
+ ,10510
+ ,5290
+ ,92.18
+ ,7480
+ ,2440
+ ,10230
+ ,5310
+ ,94.99
+ ,10450
+ ,2620
+ ,10940
+ ,4790
+ ,103.64
+ ,6360
+ ,2270
+ ,9230
+ ,3630
+ ,109.07
+ ,6710
+ ,2540
+ ,10320
+ ,2820
+ ,122.8
+ ,6200
+ ,3060
+ ,9590
+ ,2770
+ ,132.32
+ ,4490
+ ,3730
+ ,9980
+ ,2820
+ ,132.72
+ ,3480
+ ,3450
+ ,9630
+ ,3750
+ ,113.24
+ ,2520
+ ,3220
+ ,9520
+ ,3100
+ ,97.23
+ ,1920
+ ,2980
+ ,9150
+ ,4350
+ ,71.58
+ ,2010
+ ,2470
+ ,9490
+ ,4050
+ ,52.45
+ ,1950
+ ,2240
+ ,10090
+ ,6000
+ ,39.95
+ ,2240
+ ,1970
+ ,9570
+ ,10630
+ ,43.44
+ ,2370
+ ,1860
+ ,8870
+ ,3750
+ ,43.32
+ ,2840
+ ,2200
+ ,8270
+ ,3840
+ ,46.54
+ ,2700
+ ,2000
+ ,7530
+ ,4200
+ ,50.18
+ ,2980
+ ,1590
+ ,10240
+ ,2610
+ ,57.3
+ ,3290
+ ,2280
+ ,10590
+ ,2610
+ ,68.61
+ ,3300
+ ,2830
+ ,9440
+ ,2530
+ ,64.44
+ ,3000
+ ,3060
+ ,10620
+ ,3090
+ ,72.51
+ ,2330
+ ,3320
+ ,11470
+ ,4310
+ ,67.65
+ ,2190
+ ,2680
+ ,10680
+ ,4190
+ ,72.77
+ ,1970
+ ,2470
+ ,11130
+ ,3790
+ ,76.66
+ ,2170
+ ,2500
+ ,12390
+ ,3910
+ ,74.46
+ ,2830
+ ,2170
+ ,10920
+ ,4890
+ ,76.17
+ ,3190
+ ,2070
+ ,10320
+ ,4970
+ ,73.75
+ ,3550
+ ,2380
+ ,10810
+ ,5550
+ ,78.83
+ ,3240
+ ,2480
+ ,10280
+ ,4730
+ ,84.82
+ ,3450
+ ,2350
+ ,11790
+ ,4580
+ ,75.95
+ ,3570
+ ,2610
+ ,13290
+ ,2500
+ ,74.76
+ ,3230
+ ,3410
+ ,11740
+ ,2630
+ ,75.58
+ ,3260
+ ,3380
+ ,11320
+ ,4300
+ ,77.04
+ ,2700
+ ,2720
+ ,11930
+ ,3750
+ ,77.84)
+ ,dim=c(5
+ ,69)
+ ,dimnames=list(c('Garnalen'
+ ,'Kabeljauw'
+ ,'Tong'
+ ,'Zeeduivel'
+ ,'Olie')
+ ,1:69))
> y <- array(NA,dim=c(5,69),dimnames=list(c('Garnalen','Kabeljauw','Tong','Zeeduivel','Olie'),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 = 'Do not include Seasonal 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 Kabeljauw Tong Zeeduivel Olie t
1 3010 2590 11290 4700 44.51 1
2 2910 2080 11620 4960 45.48 2
3 3840 2640 10790 4880 53.10 3
4 3580 3000 8380 4090 51.88 4
5 3140 2350 9370 3450 48.65 5
6 3550 2220 10090 3020 54.35 6
7 3250 2540 11130 3070 57.52 7
8 2820 2700 10530 3720 63.98 8
9 2260 2580 10490 3750 62.91 9
10 2060 2420 10570 3910 58.54 10
11 2120 2090 11170 4120 55.24 11
12 2210 2000 11610 4780 56.86 12
13 2190 1860 10920 3070 62.99 13
14 2180 1980 11570 4100 60.21 14
15 2350 2690 12960 3900 62.06 15
16 2440 3040 11190 3020 70.26 16
17 2370 2450 11920 3220 69.78 17
18 2440 2650 14930 4030 68.56 18
19 2610 2710 14520 4210 73.67 19
20 3040 3230 12970 4510 73.23 20
21 3190 3160 13870 4320 61.96 21
22 3120 3040 13250 3890 57.81 22
23 3170 2630 12760 7280 58.76 23
24 3600 2730 14050 9640 62.47 24
25 3420 2830 14660 5680 53.68 25
26 3650 2320 15010 6320 57.56 26
27 4180 2410 15020 5820 62.05 27
28 2960 3080 13090 4890 67.49 28
29 2710 2260 13190 3320 67.21 29
30 2950 2300 11390 2930 71.05 30
31 3030 3600 10110 3530 76.93 31
32 3770 3380 8240 3690 70.76 32
33 4740 3670 7920 3750 77.17 33
34 4450 3040 7700 3330 82.34 34
35 5550 2840 7920 4790 92.41 35
36 5580 2810 8130 5990 90.93 36
37 5890 2980 10510 5290 92.18 37
38 7480 2440 10230 5310 94.99 38
39 10450 2620 10940 4790 103.64 39
40 6360 2270 9230 3630 109.07 40
41 6710 2540 10320 2820 122.80 41
42 6200 3060 9590 2770 132.32 42
43 4490 3730 9980 2820 132.72 43
44 3480 3450 9630 3750 113.24 44
45 2520 3220 9520 3100 97.23 45
46 1920 2980 9150 4350 71.58 46
47 2010 2470 9490 4050 52.45 47
48 1950 2240 10090 6000 39.95 48
49 2240 1970 9570 10630 43.44 49
50 2370 1860 8870 3750 43.32 50
51 2840 2200 8270 3840 46.54 51
52 2700 2000 7530 4200 50.18 52
53 2980 1590 10240 2610 57.30 53
54 3290 2280 10590 2610 68.61 54
55 3300 2830 9440 2530 64.44 55
56 3000 3060 10620 3090 72.51 56
57 2330 3320 11470 4310 67.65 57
58 2190 2680 10680 4190 72.77 58
59 1970 2470 11130 3790 76.66 59
60 2170 2500 12390 3910 74.46 60
61 2830 2170 10920 4890 76.17 61
62 3190 2070 10320 4970 73.75 62
63 3550 2380 10810 5550 78.83 63
64 3240 2480 10280 4730 84.82 64
65 3450 2350 11790 4580 75.95 65
66 3570 2610 13290 2500 74.76 66
67 3230 3410 11740 2630 75.58 67
68 3260 3380 11320 4300 77.04 68
69 2700 2720 11930 3750 77.84 69
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Kabeljauw Tong Zeeduivel Olie t
1054.7386 -0.3900 -0.1163 0.2447 58.4448 -15.5935
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1959.2 -738.7 -112.3 618.3 5068.2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1054.73860 1242.57471 0.849 0.3992
Kabeljauw -0.38995 0.31373 -1.243 0.2185
Tong -0.11631 0.07456 -1.560 0.1238
Zeeduivel 0.24471 0.09764 2.506 0.0148 *
Olie 58.44479 8.11064 7.206 8.74e-10 ***
t -15.59350 7.10635 -2.194 0.0319 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1091 on 63 degrees of freedom
Multiple R-squared: 0.4944, Adjusted R-squared: 0.4543
F-statistic: 12.32 on 5 and 63 DF, p-value: 2.404e-08
> 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,] 3.314870e-02 6.629740e-02 9.668513e-01
[2,] 2.477426e-02 4.954852e-02 9.752257e-01
[3,] 1.873357e-02 3.746715e-02 9.812664e-01
[4,] 1.252437e-02 2.504874e-02 9.874756e-01
[5,] 4.666444e-03 9.332888e-03 9.953336e-01
[6,] 2.585220e-03 5.170440e-03 9.974148e-01
[7,] 1.374056e-03 2.748112e-03 9.986259e-01
[8,] 4.812781e-04 9.625561e-04 9.995187e-01
[9,] 2.392548e-04 4.785096e-04 9.997607e-01
[10,] 1.198284e-04 2.396569e-04 9.998802e-01
[11,] 8.165210e-05 1.633042e-04 9.999183e-01
[12,] 8.213715e-05 1.642743e-04 9.999179e-01
[13,] 4.495473e-05 8.990947e-05 9.999550e-01
[14,] 1.606783e-05 3.213566e-05 9.999839e-01
[15,] 1.249465e-05 2.498931e-05 9.999875e-01
[16,] 6.002023e-06 1.200405e-05 9.999940e-01
[17,] 2.265912e-06 4.531824e-06 9.999977e-01
[18,] 3.316413e-06 6.632826e-06 9.999967e-01
[19,] 2.092304e-05 4.184608e-05 9.999791e-01
[20,] 8.638199e-06 1.727640e-05 9.999914e-01
[21,] 5.375811e-06 1.075162e-05 9.999946e-01
[22,] 8.340663e-06 1.668133e-05 9.999917e-01
[23,] 5.470996e-06 1.094199e-05 9.999945e-01
[24,] 2.613062e-06 5.226124e-06 9.999974e-01
[25,] 7.913492e-06 1.582698e-05 9.999921e-01
[26,] 1.185685e-05 2.371370e-05 9.999881e-01
[27,] 9.354490e-05 1.870898e-04 9.999065e-01
[28,] 8.398885e-05 1.679777e-04 9.999160e-01
[29,] 1.956964e-04 3.913928e-04 9.998043e-01
[30,] 4.553591e-03 9.107183e-03 9.954464e-01
[31,] 9.685187e-01 6.296256e-02 3.148128e-02
[32,] 9.771529e-01 4.569411e-02 2.284706e-02
[33,] 9.938726e-01 1.225471e-02 6.127357e-03
[34,] 9.996206e-01 7.587220e-04 3.793610e-04
[35,] 9.999099e-01 1.801527e-04 9.007634e-05
[36,] 9.999772e-01 4.550801e-05 2.275400e-05
[37,] 9.999785e-01 4.291413e-05 2.145706e-05
[38,] 9.999650e-01 7.009466e-05 3.504733e-05
[39,] 9.999151e-01 1.697437e-04 8.487184e-05
[40,] 9.997962e-01 4.075173e-04 2.037586e-04
[41,] 9.998125e-01 3.749955e-04 1.874978e-04
[42,] 9.995644e-01 8.711163e-04 4.355581e-04
[43,] 9.989403e-01 2.119340e-03 1.059670e-03
[44,] 9.980637e-01 3.872684e-03 1.936342e-03
[45,] 9.960380e-01 7.923936e-03 3.961968e-03
[46,] 9.929322e-01 1.413567e-02 7.067833e-03
[47,] 9.857544e-01 2.849120e-02 1.424560e-02
[48,] 9.886454e-01 2.270916e-02 1.135458e-02
[49,] 9.741006e-01 5.179879e-02 2.589939e-02
[50,] 9.406518e-01 1.186964e-01 5.934818e-02
[51,] 9.067891e-01 1.864217e-01 9.321087e-02
[52,] 9.432040e-01 1.135919e-01 5.679597e-02
> postscript(file="/var/www/html/rcomp/tmp/1n2qy1293212701.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/2gb711293212701.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/3gb711293212701.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/4gb711293212701.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/5gb711293212701.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
542.47344 177.25685 818.91509 699.20704 481.86628 712.59776
7 8 9 10 11 12
476.43373 -481.97924 -1022.63802 -1043.88197 -885.70839 -1020.22203
13 14 15 16 17 18
-1099.29263 -1060.87542 -495.92458 -723.62008 -944.08168 -557.31487
19 20 21 22 23 24
-738.71141 -318.31837 629.82445 804.28137 -232.07974 -391.79107
25 26 27 28 29 30
1036.51874 740.56512 1182.35433 -75.62222 -217.60519 -284.76339
31 32 33 34 35 36
-321.58531 452.16973 1124.31737 379.26879 496.64738 347.81650
37 38 39 40 41 42
1114.75924 2308.08637 5068.15227 624.87722 618.30249 -302.39280
43 44 45 46 47 48
-1725.78203 -1959.15788 -1911.28703 -1439.09295 -301.36989 -112.29996
49 50 51 52 53 54
-1309.44389 402.43925 740.61551 151.31493 575.18501 549.54519
55 56 57 58 59 60
919.14880 252.99188 -215.66453 -951.39935 -1314.82409 -841.76782
61 62 63 64 65 66
-805.58902 -416.91717 -302.27577 -768.75513 136.88366 1126.87137
67 68 69
854.41160 345.46521 -297.52905
> postscript(file="/var/www/html/rcomp/tmp/6qkp41293212701.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 542.47344 NA
1 177.25685 542.47344
2 818.91509 177.25685
3 699.20704 818.91509
4 481.86628 699.20704
5 712.59776 481.86628
6 476.43373 712.59776
7 -481.97924 476.43373
8 -1022.63802 -481.97924
9 -1043.88197 -1022.63802
10 -885.70839 -1043.88197
11 -1020.22203 -885.70839
12 -1099.29263 -1020.22203
13 -1060.87542 -1099.29263
14 -495.92458 -1060.87542
15 -723.62008 -495.92458
16 -944.08168 -723.62008
17 -557.31487 -944.08168
18 -738.71141 -557.31487
19 -318.31837 -738.71141
20 629.82445 -318.31837
21 804.28137 629.82445
22 -232.07974 804.28137
23 -391.79107 -232.07974
24 1036.51874 -391.79107
25 740.56512 1036.51874
26 1182.35433 740.56512
27 -75.62222 1182.35433
28 -217.60519 -75.62222
29 -284.76339 -217.60519
30 -321.58531 -284.76339
31 452.16973 -321.58531
32 1124.31737 452.16973
33 379.26879 1124.31737
34 496.64738 379.26879
35 347.81650 496.64738
36 1114.75924 347.81650
37 2308.08637 1114.75924
38 5068.15227 2308.08637
39 624.87722 5068.15227
40 618.30249 624.87722
41 -302.39280 618.30249
42 -1725.78203 -302.39280
43 -1959.15788 -1725.78203
44 -1911.28703 -1959.15788
45 -1439.09295 -1911.28703
46 -301.36989 -1439.09295
47 -112.29996 -301.36989
48 -1309.44389 -112.29996
49 402.43925 -1309.44389
50 740.61551 402.43925
51 151.31493 740.61551
52 575.18501 151.31493
53 549.54519 575.18501
54 919.14880 549.54519
55 252.99188 919.14880
56 -215.66453 252.99188
57 -951.39935 -215.66453
58 -1314.82409 -951.39935
59 -841.76782 -1314.82409
60 -805.58902 -841.76782
61 -416.91717 -805.58902
62 -302.27577 -416.91717
63 -768.75513 -302.27577
64 136.88366 -768.75513
65 1126.87137 136.88366
66 854.41160 1126.87137
67 345.46521 854.41160
68 -297.52905 345.46521
69 NA -297.52905
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 177.25685 542.47344
[2,] 818.91509 177.25685
[3,] 699.20704 818.91509
[4,] 481.86628 699.20704
[5,] 712.59776 481.86628
[6,] 476.43373 712.59776
[7,] -481.97924 476.43373
[8,] -1022.63802 -481.97924
[9,] -1043.88197 -1022.63802
[10,] -885.70839 -1043.88197
[11,] -1020.22203 -885.70839
[12,] -1099.29263 -1020.22203
[13,] -1060.87542 -1099.29263
[14,] -495.92458 -1060.87542
[15,] -723.62008 -495.92458
[16,] -944.08168 -723.62008
[17,] -557.31487 -944.08168
[18,] -738.71141 -557.31487
[19,] -318.31837 -738.71141
[20,] 629.82445 -318.31837
[21,] 804.28137 629.82445
[22,] -232.07974 804.28137
[23,] -391.79107 -232.07974
[24,] 1036.51874 -391.79107
[25,] 740.56512 1036.51874
[26,] 1182.35433 740.56512
[27,] -75.62222 1182.35433
[28,] -217.60519 -75.62222
[29,] -284.76339 -217.60519
[30,] -321.58531 -284.76339
[31,] 452.16973 -321.58531
[32,] 1124.31737 452.16973
[33,] 379.26879 1124.31737
[34,] 496.64738 379.26879
[35,] 347.81650 496.64738
[36,] 1114.75924 347.81650
[37,] 2308.08637 1114.75924
[38,] 5068.15227 2308.08637
[39,] 624.87722 5068.15227
[40,] 618.30249 624.87722
[41,] -302.39280 618.30249
[42,] -1725.78203 -302.39280
[43,] -1959.15788 -1725.78203
[44,] -1911.28703 -1959.15788
[45,] -1439.09295 -1911.28703
[46,] -301.36989 -1439.09295
[47,] -112.29996 -301.36989
[48,] -1309.44389 -112.29996
[49,] 402.43925 -1309.44389
[50,] 740.61551 402.43925
[51,] 151.31493 740.61551
[52,] 575.18501 151.31493
[53,] 549.54519 575.18501
[54,] 919.14880 549.54519
[55,] 252.99188 919.14880
[56,] -215.66453 252.99188
[57,] -951.39935 -215.66453
[58,] -1314.82409 -951.39935
[59,] -841.76782 -1314.82409
[60,] -805.58902 -841.76782
[61,] -416.91717 -805.58902
[62,] -302.27577 -416.91717
[63,] -768.75513 -302.27577
[64,] 136.88366 -768.75513
[65,] 1126.87137 136.88366
[66,] 854.41160 1126.87137
[67,] 345.46521 854.41160
[68,] -297.52905 345.46521
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 177.25685 542.47344
2 818.91509 177.25685
3 699.20704 818.91509
4 481.86628 699.20704
5 712.59776 481.86628
6 476.43373 712.59776
7 -481.97924 476.43373
8 -1022.63802 -481.97924
9 -1043.88197 -1022.63802
10 -885.70839 -1043.88197
11 -1020.22203 -885.70839
12 -1099.29263 -1020.22203
13 -1060.87542 -1099.29263
14 -495.92458 -1060.87542
15 -723.62008 -495.92458
16 -944.08168 -723.62008
17 -557.31487 -944.08168
18 -738.71141 -557.31487
19 -318.31837 -738.71141
20 629.82445 -318.31837
21 804.28137 629.82445
22 -232.07974 804.28137
23 -391.79107 -232.07974
24 1036.51874 -391.79107
25 740.56512 1036.51874
26 1182.35433 740.56512
27 -75.62222 1182.35433
28 -217.60519 -75.62222
29 -284.76339 -217.60519
30 -321.58531 -284.76339
31 452.16973 -321.58531
32 1124.31737 452.16973
33 379.26879 1124.31737
34 496.64738 379.26879
35 347.81650 496.64738
36 1114.75924 347.81650
37 2308.08637 1114.75924
38 5068.15227 2308.08637
39 624.87722 5068.15227
40 618.30249 624.87722
41 -302.39280 618.30249
42 -1725.78203 -302.39280
43 -1959.15788 -1725.78203
44 -1911.28703 -1959.15788
45 -1439.09295 -1911.28703
46 -301.36989 -1439.09295
47 -112.29996 -301.36989
48 -1309.44389 -112.29996
49 402.43925 -1309.44389
50 740.61551 402.43925
51 151.31493 740.61551
52 575.18501 151.31493
53 549.54519 575.18501
54 919.14880 549.54519
55 252.99188 919.14880
56 -215.66453 252.99188
57 -951.39935 -215.66453
58 -1314.82409 -951.39935
59 -841.76782 -1314.82409
60 -805.58902 -841.76782
61 -416.91717 -805.58902
62 -302.27577 -416.91717
63 -768.75513 -302.27577
64 136.88366 -768.75513
65 1126.87137 136.88366
66 854.41160 1126.87137
67 345.46521 854.41160
68 -297.52905 345.46521
> 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/7ulns1293212701.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/8ulns1293212701.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/9mcmd1293212701.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/10f3mx1293212701.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/11td161293212701.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/12mmjr1293212701.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/13oqd01293212702.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/14hzu31293212702.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/155irf1293212702.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/16r1ql1293212702.tab")
+ }
>
> try(system("convert tmp/1n2qy1293212701.ps tmp/1n2qy1293212701.png",intern=TRUE))
character(0)
> try(system("convert tmp/2gb711293212701.ps tmp/2gb711293212701.png",intern=TRUE))
character(0)
> try(system("convert tmp/3gb711293212701.ps tmp/3gb711293212701.png",intern=TRUE))
character(0)
> try(system("convert tmp/4gb711293212701.ps tmp/4gb711293212701.png",intern=TRUE))
character(0)
> try(system("convert tmp/5gb711293212701.ps tmp/5gb711293212701.png",intern=TRUE))
character(0)
> try(system("convert tmp/6qkp41293212701.ps tmp/6qkp41293212701.png",intern=TRUE))
character(0)
> try(system("convert tmp/7ulns1293212701.ps tmp/7ulns1293212701.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ulns1293212701.ps tmp/8ulns1293212701.png",intern=TRUE))
character(0)
> try(system("convert tmp/9mcmd1293212701.ps tmp/9mcmd1293212701.png",intern=TRUE))
character(0)
> try(system("convert tmp/10f3mx1293212701.ps tmp/10f3mx1293212701.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.744 1.750 7.394