R version 2.12.0 (2010-10-15)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(7378.00
+ ,7599.00
+ ,7869.78
+ ,5490.00
+ ,5666.00
+ ,5892.00
+ ,4082.00
+ ,4193.00
+ ,4358.00
+ ,4340.00
+ ,4463.00
+ ,4634.00
+ ,3667.00
+ ,3769.00
+ ,3906.00
+ ,3950.00
+ ,4029.00
+ ,4152.00
+ ,3753.00
+ ,3865.00
+ ,4008.00
+ ,4037.00
+ ,4140.00
+ ,4277.00
+ ,3912.00
+ ,4011.00
+ ,4149.00
+ ,3207.00
+ ,3261.00
+ ,3366.00
+ ,3892.00
+ ,3996.00
+ ,4137.00
+ ,3427.00
+ ,3518.00
+ ,3657.00
+ ,3162.00
+ ,3240.00
+ ,3347.00
+ ,2974.00
+ ,3045.00
+ ,3143.00
+ ,3032.00
+ ,3108.00
+ ,3214.53
+ ,3469.00
+ ,3563.00
+ ,3697.00
+ ,3191.00
+ ,3293.00
+ ,3409.75
+ ,2836.00
+ ,2907.00
+ ,3000.00
+ ,3076.00
+ ,3200.00
+ ,3337.00
+ ,3178.00
+ ,3250.00
+ ,3357.00
+ ,2984.00
+ ,3068.00
+ ,3178.00
+ ,2595.00
+ ,2664.00
+ ,2763.00
+ ,2764.00
+ ,2851.00
+ ,2956.00
+ ,2619.00
+ ,2673.00
+ ,2759.00
+ ,2094.00
+ ,2155.05
+ ,2240.00
+ ,2649.00
+ ,2705.00
+ ,2783.00
+ ,2303.00
+ ,2360.00
+ ,2438.00
+ ,2208.00
+ ,2263.00
+ ,2336.00
+ ,2957.00
+ ,3025.00
+ ,3124.00
+ ,1848.00
+ ,1899.00
+ ,1975.00
+ ,2468.00
+ ,2527.00
+ ,2607.00
+ ,2133.00
+ ,2172.00
+ ,2236.00
+ ,2537.00
+ ,2588.00
+ ,2669.00
+ ,2341.00
+ ,2402.00
+ ,2487.00
+ ,2334.00
+ ,2375.00
+ ,2449.00
+ ,2285.00
+ ,2341.00
+ ,2420.00
+ ,2366.00
+ ,2455.00
+ ,2551.00
+ ,2450.00
+ ,2509.00
+ ,2590.00
+ ,2484.00
+ ,2563.00
+ ,2667.00
+ ,2458.00
+ ,2537.00
+ ,2629.00
+ ,2425.00
+ ,2490.00
+ ,2582.00
+ ,2049.00
+ ,2105.00
+ ,2191.00
+ ,2037.00
+ ,2096.00
+ ,2180.00
+ ,2455.00
+ ,2548.00
+ ,2657.00
+ ,2118.00
+ ,2180.00
+ ,2267.00
+ ,2093.00
+ ,2156.00
+ ,2243.00
+ ,2097.00
+ ,2128.00
+ ,2193.00
+ ,1988.00
+ ,2041.00
+ ,2126.00
+ ,2510.00
+ ,2561.00
+ ,2641.00
+ ,2474.00
+ ,2509.00
+ ,2590.00
+ ,2202.00
+ ,2263.00
+ ,2356.00
+ ,2407.00
+ ,2464.00
+ ,2551.00
+ ,3166.00
+ ,3238.00
+ ,3334.00
+ ,2522.00
+ ,2573.00
+ ,2647.00
+ ,2486.00
+ ,2553.00
+ ,2649.00
+ ,2753.00
+ ,2817.00
+ ,2903.00
+ ,2500.00
+ ,2574.00
+ ,2668.00
+ ,2081.00
+ ,2142.00
+ ,2223.00
+ ,2525.00
+ ,2592.00
+ ,2685.00
+ ,2165.00
+ ,2245.00
+ ,2334.00
+ ,2299.00
+ ,2346.00
+ ,2409.00
+ ,2361.00
+ ,2433.00
+ ,2532.00
+ ,2613.00
+ ,2672.00
+ ,2757.00
+ ,2643.00
+ ,2705.00
+ ,2785.00
+ ,2248.00
+ ,2316.00
+ ,2412.00
+ ,2389.00
+ ,2454.00
+ ,2549.00
+ ,3135.00
+ ,3199.00
+ ,3303.00
+ ,2159.00
+ ,2230.00
+ ,2328.00
+ ,1897.00
+ ,1965.00
+ ,2047.00
+ ,2022.00
+ ,2090.00
+ ,2175.00
+ ,2106.00
+ ,2162.00
+ ,2249.00
+ ,2028.00
+ ,2078.00
+ ,2149.00
+ ,2212.00
+ ,2274.00
+ ,2373.00
+ ,2169.00
+ ,2247.00
+ ,2332.00
+ ,2252.00
+ ,2306.00
+ ,2370.00)
+ ,dim=c(3
+ ,75)
+ ,dimnames=list(c('Loon2006'
+ ,'Loon2007'
+ ,'Loon2008')
+ ,1:75))
> y <- array(NA,dim=c(3,75),dimnames=list(c('Loon2006','Loon2007','Loon2008'),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 = 'Do not include Seasonal Dummies'
> par1 = '3'
> #'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
> 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
Loon2008 Loon2006 Loon2007
1 7869.78 7378 7599.00
2 5892.00 5490 5666.00
3 4358.00 4082 4193.00
4 4634.00 4340 4463.00
5 3906.00 3667 3769.00
6 4152.00 3950 4029.00
7 4008.00 3753 3865.00
8 4277.00 4037 4140.00
9 4149.00 3912 4011.00
10 3366.00 3207 3261.00
11 4137.00 3892 3996.00
12 3657.00 3427 3518.00
13 3347.00 3162 3240.00
14 3143.00 2974 3045.00
15 3214.53 3032 3108.00
16 3697.00 3469 3563.00
17 3409.75 3191 3293.00
18 3000.00 2836 2907.00
19 3337.00 3076 3200.00
20 3357.00 3178 3250.00
21 3178.00 2984 3068.00
22 2763.00 2595 2664.00
23 2956.00 2764 2851.00
24 2759.00 2619 2673.00
25 2240.00 2094 2155.05
26 2783.00 2649 2705.00
27 2438.00 2303 2360.00
28 2336.00 2208 2263.00
29 3124.00 2957 3025.00
30 1975.00 1848 1899.00
31 2607.00 2468 2527.00
32 2236.00 2133 2172.00
33 2669.00 2537 2588.00
34 2487.00 2341 2402.00
35 2449.00 2334 2375.00
36 2420.00 2285 2341.00
37 2551.00 2366 2455.00
38 2590.00 2450 2509.00
39 2667.00 2484 2563.00
40 2629.00 2458 2537.00
41 2582.00 2425 2490.00
42 2191.00 2049 2105.00
43 2180.00 2037 2096.00
44 2657.00 2455 2548.00
45 2267.00 2118 2180.00
46 2243.00 2093 2156.00
47 2193.00 2097 2128.00
48 2126.00 1988 2041.00
49 2641.00 2510 2561.00
50 2590.00 2474 2509.00
51 2356.00 2202 2263.00
52 2551.00 2407 2464.00
53 3334.00 3166 3238.00
54 2647.00 2522 2573.00
55 2649.00 2486 2553.00
56 2903.00 2753 2817.00
57 2668.00 2500 2574.00
58 2223.00 2081 2142.00
59 2685.00 2525 2592.00
60 2334.00 2165 2245.00
61 2409.00 2299 2346.00
62 2532.00 2361 2433.00
63 2757.00 2613 2672.00
64 2785.00 2643 2705.00
65 2412.00 2248 2316.00
66 2549.00 2389 2454.00
67 3303.00 3135 3199.00
68 2328.00 2159 2230.00
69 2047.00 1897 1965.00
70 2175.00 2022 2090.00
71 2249.00 2106 2162.00
72 2149.00 2028 2078.00
73 2373.00 2212 2274.00
74 2332.00 2169 2247.00
75 2370.00 2252 2306.00
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Loon2006 Loon2007
5.016 -0.661 1.678
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-15.3989 -5.0722 -0.2747 4.2976 16.2642
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.01557 2.76061 1.817 0.0734 .
Loon2006 -0.66103 0.06037 -10.950 <2e-16 ***
Loon2007 1.67780 0.05865 28.605 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.038 on 72 degrees of freedom
Multiple R-squared: 0.9999, Adjusted R-squared: 0.9999
F-statistic: 6.337e+05 on 2 and 72 DF, p-value: < 2.2e-16
> 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.725908800 0.54818240 0.27409120
[2,] 0.837191750 0.32561650 0.16280825
[3,] 0.838942591 0.32211482 0.16105741
[4,] 0.758572668 0.48285466 0.24142733
[5,] 0.795509604 0.40898079 0.20449040
[6,] 0.740094584 0.51981083 0.25990542
[7,] 0.835546423 0.32890715 0.16445358
[8,] 0.870526625 0.25894675 0.12947338
[9,] 0.891815590 0.21636882 0.10818441
[10,] 0.864284694 0.27143061 0.13571531
[11,] 0.873855224 0.25228955 0.12614478
[12,] 0.940077654 0.11984469 0.05992235
[13,] 0.944295691 0.11140862 0.05570431
[14,] 0.920686721 0.15862656 0.07931328
[15,] 0.897799981 0.20440004 0.10220002
[16,] 0.865847331 0.26830534 0.13415267
[17,] 0.841103423 0.31779315 0.15889658
[18,] 0.805526454 0.38894709 0.19447355
[19,] 0.753953108 0.49209378 0.24604689
[20,] 0.708714237 0.58257153 0.29128576
[21,] 0.751881482 0.49623704 0.24811852
[22,] 0.712583531 0.57483294 0.28741647
[23,] 0.701206390 0.59758722 0.29879361
[24,] 0.647178619 0.70564276 0.35282138
[25,] 0.622650912 0.75469818 0.37734909
[26,] 0.600357734 0.79928453 0.39964227
[27,] 0.564736903 0.87052619 0.43526310
[28,] 0.496131388 0.99226278 0.50386861
[29,] 0.428082151 0.85616430 0.57191785
[30,] 0.366269680 0.73253936 0.63373032
[31,] 0.311656524 0.62331305 0.68834348
[32,] 0.308219529 0.61643906 0.69178047
[33,] 0.273292048 0.54658410 0.72670795
[34,] 0.257736960 0.51547392 0.74226304
[35,] 0.242491855 0.48498371 0.75750815
[36,] 0.202952855 0.40590571 0.79704715
[37,] 0.226313484 0.45262697 0.77368652
[38,] 0.198427769 0.39685554 0.80157223
[39,] 0.156852710 0.31370542 0.84314729
[40,] 0.131455223 0.26291045 0.86854478
[41,] 0.106896732 0.21379346 0.89310327
[42,] 0.080791858 0.16158372 0.91920814
[43,] 0.100530296 0.20106059 0.89946970
[44,] 0.073983174 0.14796635 0.92601683
[45,] 0.118116303 0.23623261 0.88188370
[46,] 0.151288011 0.30257602 0.84871199
[47,] 0.127021163 0.25404233 0.87297884
[48,] 0.149837183 0.29967437 0.85016282
[49,] 0.135784571 0.27156914 0.86421543
[50,] 0.110786931 0.22157386 0.88921307
[51,] 0.109455284 0.21891057 0.89054472
[52,] 0.083221967 0.16644393 0.91677803
[53,] 0.056384787 0.11276957 0.94361521
[54,] 0.036617625 0.07323525 0.96338237
[55,] 0.045075123 0.09015025 0.95492488
[56,] 0.058675658 0.11735132 0.94132434
[57,] 0.042957120 0.08591424 0.95704288
[58,] 0.028050037 0.05610007 0.97194996
[59,] 0.054241564 0.10848313 0.94575844
[60,] 0.041642056 0.08328411 0.95835794
[61,] 0.027599193 0.05519839 0.97240081
[62,] 0.014486328 0.02897266 0.98551367
[63,] 0.015167681 0.03033536 0.98483232
[64,] 0.006653031 0.01330606 0.99334697
> postscript(file="/var/www/rcomp/tmp/1qjfp1321714989.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/rcomp/tmp/2cldi1321714989.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/rcomp/tmp/3xmxd1321714989.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/rcomp/tmp/4u8km1321714989.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/rcomp/tmp/5eoqx1321714989.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 75
Frequency = 1
1 2 3 4 5 6
-7.81141046 9.58486071 16.26422634 9.80209590 1.32694254 -1.83147017
7 8 9 10 11 12
-0.89392565 -5.55837309 0.24997565 9.57898702 0.19649622 14.80928181
13 14 15 16 17 18
-3.93334423 -5.03462846 -0.86671843 7.07123301 -10.93717777 -7.71939443
19 20 21 22 23 24
-3.66950671 -0.13495244 -2.01385259 3.67948482 -5.35631215 0.44389721
25 26 27 28 29 30
3.42333370 -9.41502277 -4.28796477 -6.33853600 -1.71600611 5.41242141
31 32 33 34 35 36
-6.41178760 -3.23539073 -1.14697166 -0.63670945 2.03681234 -2.30817360
37 38 39 40 41 42
-9.03465125 -5.10980055 3.76369708 -7.80009940 2.24280118 8.65121212
43 44 45 46 47 48
4.81912446 -0.23902590 4.42677040 4.16839273 3.79101631 10.70802503
49 50 51 52 53 54
-1.69399109 10.75484891 9.69530163 2.96722118 -10.93362774 -7.89531579
55 56 57 58 59 60
3.86379240 -8.58226977 -3.11571525 -0.27467433 0.20948713 -6.56222548
61 62 63 64 65 66
-12.44281534 5.57190407 -3.84446103 -11.38118514 7.17892812 5.84677528
67 68 69 70 71 72
3.00889401 8.63867394 -0.93232447 -0.02945674 8.69491981 -1.92964494
73 74 75
14.84972693 -9.27372548 -15.39892244
> postscript(file="/var/www/rcomp/tmp/6jeuk1321714989.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 -7.81141046 NA
1 9.58486071 -7.81141046
2 16.26422634 9.58486071
3 9.80209590 16.26422634
4 1.32694254 9.80209590
5 -1.83147017 1.32694254
6 -0.89392565 -1.83147017
7 -5.55837309 -0.89392565
8 0.24997565 -5.55837309
9 9.57898702 0.24997565
10 0.19649622 9.57898702
11 14.80928181 0.19649622
12 -3.93334423 14.80928181
13 -5.03462846 -3.93334423
14 -0.86671843 -5.03462846
15 7.07123301 -0.86671843
16 -10.93717777 7.07123301
17 -7.71939443 -10.93717777
18 -3.66950671 -7.71939443
19 -0.13495244 -3.66950671
20 -2.01385259 -0.13495244
21 3.67948482 -2.01385259
22 -5.35631215 3.67948482
23 0.44389721 -5.35631215
24 3.42333370 0.44389721
25 -9.41502277 3.42333370
26 -4.28796477 -9.41502277
27 -6.33853600 -4.28796477
28 -1.71600611 -6.33853600
29 5.41242141 -1.71600611
30 -6.41178760 5.41242141
31 -3.23539073 -6.41178760
32 -1.14697166 -3.23539073
33 -0.63670945 -1.14697166
34 2.03681234 -0.63670945
35 -2.30817360 2.03681234
36 -9.03465125 -2.30817360
37 -5.10980055 -9.03465125
38 3.76369708 -5.10980055
39 -7.80009940 3.76369708
40 2.24280118 -7.80009940
41 8.65121212 2.24280118
42 4.81912446 8.65121212
43 -0.23902590 4.81912446
44 4.42677040 -0.23902590
45 4.16839273 4.42677040
46 3.79101631 4.16839273
47 10.70802503 3.79101631
48 -1.69399109 10.70802503
49 10.75484891 -1.69399109
50 9.69530163 10.75484891
51 2.96722118 9.69530163
52 -10.93362774 2.96722118
53 -7.89531579 -10.93362774
54 3.86379240 -7.89531579
55 -8.58226977 3.86379240
56 -3.11571525 -8.58226977
57 -0.27467433 -3.11571525
58 0.20948713 -0.27467433
59 -6.56222548 0.20948713
60 -12.44281534 -6.56222548
61 5.57190407 -12.44281534
62 -3.84446103 5.57190407
63 -11.38118514 -3.84446103
64 7.17892812 -11.38118514
65 5.84677528 7.17892812
66 3.00889401 5.84677528
67 8.63867394 3.00889401
68 -0.93232447 8.63867394
69 -0.02945674 -0.93232447
70 8.69491981 -0.02945674
71 -1.92964494 8.69491981
72 14.84972693 -1.92964494
73 -9.27372548 14.84972693
74 -15.39892244 -9.27372548
75 NA -15.39892244
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 9.58486071 -7.81141046
[2,] 16.26422634 9.58486071
[3,] 9.80209590 16.26422634
[4,] 1.32694254 9.80209590
[5,] -1.83147017 1.32694254
[6,] -0.89392565 -1.83147017
[7,] -5.55837309 -0.89392565
[8,] 0.24997565 -5.55837309
[9,] 9.57898702 0.24997565
[10,] 0.19649622 9.57898702
[11,] 14.80928181 0.19649622
[12,] -3.93334423 14.80928181
[13,] -5.03462846 -3.93334423
[14,] -0.86671843 -5.03462846
[15,] 7.07123301 -0.86671843
[16,] -10.93717777 7.07123301
[17,] -7.71939443 -10.93717777
[18,] -3.66950671 -7.71939443
[19,] -0.13495244 -3.66950671
[20,] -2.01385259 -0.13495244
[21,] 3.67948482 -2.01385259
[22,] -5.35631215 3.67948482
[23,] 0.44389721 -5.35631215
[24,] 3.42333370 0.44389721
[25,] -9.41502277 3.42333370
[26,] -4.28796477 -9.41502277
[27,] -6.33853600 -4.28796477
[28,] -1.71600611 -6.33853600
[29,] 5.41242141 -1.71600611
[30,] -6.41178760 5.41242141
[31,] -3.23539073 -6.41178760
[32,] -1.14697166 -3.23539073
[33,] -0.63670945 -1.14697166
[34,] 2.03681234 -0.63670945
[35,] -2.30817360 2.03681234
[36,] -9.03465125 -2.30817360
[37,] -5.10980055 -9.03465125
[38,] 3.76369708 -5.10980055
[39,] -7.80009940 3.76369708
[40,] 2.24280118 -7.80009940
[41,] 8.65121212 2.24280118
[42,] 4.81912446 8.65121212
[43,] -0.23902590 4.81912446
[44,] 4.42677040 -0.23902590
[45,] 4.16839273 4.42677040
[46,] 3.79101631 4.16839273
[47,] 10.70802503 3.79101631
[48,] -1.69399109 10.70802503
[49,] 10.75484891 -1.69399109
[50,] 9.69530163 10.75484891
[51,] 2.96722118 9.69530163
[52,] -10.93362774 2.96722118
[53,] -7.89531579 -10.93362774
[54,] 3.86379240 -7.89531579
[55,] -8.58226977 3.86379240
[56,] -3.11571525 -8.58226977
[57,] -0.27467433 -3.11571525
[58,] 0.20948713 -0.27467433
[59,] -6.56222548 0.20948713
[60,] -12.44281534 -6.56222548
[61,] 5.57190407 -12.44281534
[62,] -3.84446103 5.57190407
[63,] -11.38118514 -3.84446103
[64,] 7.17892812 -11.38118514
[65,] 5.84677528 7.17892812
[66,] 3.00889401 5.84677528
[67,] 8.63867394 3.00889401
[68,] -0.93232447 8.63867394
[69,] -0.02945674 -0.93232447
[70,] 8.69491981 -0.02945674
[71,] -1.92964494 8.69491981
[72,] 14.84972693 -1.92964494
[73,] -9.27372548 14.84972693
[74,] -15.39892244 -9.27372548
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 9.58486071 -7.81141046
2 16.26422634 9.58486071
3 9.80209590 16.26422634
4 1.32694254 9.80209590
5 -1.83147017 1.32694254
6 -0.89392565 -1.83147017
7 -5.55837309 -0.89392565
8 0.24997565 -5.55837309
9 9.57898702 0.24997565
10 0.19649622 9.57898702
11 14.80928181 0.19649622
12 -3.93334423 14.80928181
13 -5.03462846 -3.93334423
14 -0.86671843 -5.03462846
15 7.07123301 -0.86671843
16 -10.93717777 7.07123301
17 -7.71939443 -10.93717777
18 -3.66950671 -7.71939443
19 -0.13495244 -3.66950671
20 -2.01385259 -0.13495244
21 3.67948482 -2.01385259
22 -5.35631215 3.67948482
23 0.44389721 -5.35631215
24 3.42333370 0.44389721
25 -9.41502277 3.42333370
26 -4.28796477 -9.41502277
27 -6.33853600 -4.28796477
28 -1.71600611 -6.33853600
29 5.41242141 -1.71600611
30 -6.41178760 5.41242141
31 -3.23539073 -6.41178760
32 -1.14697166 -3.23539073
33 -0.63670945 -1.14697166
34 2.03681234 -0.63670945
35 -2.30817360 2.03681234
36 -9.03465125 -2.30817360
37 -5.10980055 -9.03465125
38 3.76369708 -5.10980055
39 -7.80009940 3.76369708
40 2.24280118 -7.80009940
41 8.65121212 2.24280118
42 4.81912446 8.65121212
43 -0.23902590 4.81912446
44 4.42677040 -0.23902590
45 4.16839273 4.42677040
46 3.79101631 4.16839273
47 10.70802503 3.79101631
48 -1.69399109 10.70802503
49 10.75484891 -1.69399109
50 9.69530163 10.75484891
51 2.96722118 9.69530163
52 -10.93362774 2.96722118
53 -7.89531579 -10.93362774
54 3.86379240 -7.89531579
55 -8.58226977 3.86379240
56 -3.11571525 -8.58226977
57 -0.27467433 -3.11571525
58 0.20948713 -0.27467433
59 -6.56222548 0.20948713
60 -12.44281534 -6.56222548
61 5.57190407 -12.44281534
62 -3.84446103 5.57190407
63 -11.38118514 -3.84446103
64 7.17892812 -11.38118514
65 5.84677528 7.17892812
66 3.00889401 5.84677528
67 8.63867394 3.00889401
68 -0.93232447 8.63867394
69 -0.02945674 -0.93232447
70 8.69491981 -0.02945674
71 -1.92964494 8.69491981
72 14.84972693 -1.92964494
73 -9.27372548 14.84972693
74 -15.39892244 -9.27372548
> 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/rcomp/tmp/7irh91321714989.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/rcomp/tmp/8s2ni1321714989.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/rcomp/tmp/9m38z1321714989.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/rcomp/tmp/10omhm1321714989.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/118fz31321714989.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/rcomp/tmp/12okq21321714989.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/rcomp/tmp/138g0s1321714989.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/rcomp/tmp/14uucc1321714989.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/rcomp/tmp/15mtme1321714989.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/rcomp/tmp/16v9xi1321714989.tab")
+ }
>
> try(system("convert tmp/1qjfp1321714989.ps tmp/1qjfp1321714989.png",intern=TRUE))
character(0)
> try(system("convert tmp/2cldi1321714989.ps tmp/2cldi1321714989.png",intern=TRUE))
character(0)
> try(system("convert tmp/3xmxd1321714989.ps tmp/3xmxd1321714989.png",intern=TRUE))
character(0)
> try(system("convert tmp/4u8km1321714989.ps tmp/4u8km1321714989.png",intern=TRUE))
character(0)
> try(system("convert tmp/5eoqx1321714989.ps tmp/5eoqx1321714989.png",intern=TRUE))
character(0)
> try(system("convert tmp/6jeuk1321714989.ps tmp/6jeuk1321714989.png",intern=TRUE))
character(0)
> try(system("convert tmp/7irh91321714989.ps tmp/7irh91321714989.png",intern=TRUE))
character(0)
> try(system("convert tmp/8s2ni1321714989.ps tmp/8s2ni1321714989.png",intern=TRUE))
character(0)
> try(system("convert tmp/9m38z1321714989.ps tmp/9m38z1321714989.png",intern=TRUE))
character(0)
> try(system("convert tmp/10omhm1321714989.ps tmp/10omhm1321714989.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.470 0.310 5.043