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(20604.6,2.05,18714.9,2.03,18492.6,2.04,18183.6,2.03,19435.1,2.01,22686.8,2.01,20396.7,2.01,19233.6,2.01,22751,2.01,19864,2.01,17165.4,2.02,22309.7,2.02,21786.3,2.03,21927.6,2.05,20957.9,2.08,19726,2.07,21315.7,2.06,24771.5,2.05,22592.4,2.05,21942.1,2.05,23973.7,2.05,20815.7,2.05,19931.4,2.06,24436.8,2.06,22838.7,2.07,24465.3,2.07,23007.3,2.3,22720.8,2.31,23045.7,2.31,27198.5,2.53,22401.9,2.58,25122.7,2.59,26100.5,2.73,22904.9,2.82,22040.4,3,25981.5,3.04,26157.1,3.23,25975.4,3.32,22589.8,3.49,25370.4,3.57,25091.1,3.56,28760.9,3.72,24325.9,3.82,25821.7,3.82,27645.7,3.98,26296.9,4.06,24141.5,4.08,27268.1,4.19,29060.3,4.16,28226.4,4.17,23268.5,4.21,26938.2,4.21,27217.5,4.17,27540.5,4.19,29167.6,4.25,26671.5,4.25,30184,4.2,28422.3,4.33,23774.3,4.41,29601,4.56,28523.6,5.18,23622,3.42,21320.3,2.71,20423.6,2.29,21174.9,2,23050.2,1.64,21202.9,1.3,20476.4,1.08,23173.3,1,22468,1,19842.7,1),dim=c(2,71),dimnames=list(c('Y','X'),1:71))
> y <- array(NA,dim=c(2,71),dimnames=list(c('Y','X'),1:71))
> 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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 20604.6 2.05 1 0 0 0 0 0 0 0 0 0 0 1
2 18714.9 2.03 0 1 0 0 0 0 0 0 0 0 0 2
3 18492.6 2.04 0 0 1 0 0 0 0 0 0 0 0 3
4 18183.6 2.03 0 0 0 1 0 0 0 0 0 0 0 4
5 19435.1 2.01 0 0 0 0 1 0 0 0 0 0 0 5
6 22686.8 2.01 0 0 0 0 0 1 0 0 0 0 0 6
7 20396.7 2.01 0 0 0 0 0 0 1 0 0 0 0 7
8 19233.6 2.01 0 0 0 0 0 0 0 1 0 0 0 8
9 22751.0 2.01 0 0 0 0 0 0 0 0 1 0 0 9
10 19864.0 2.01 0 0 0 0 0 0 0 0 0 1 0 10
11 17165.4 2.02 0 0 0 0 0 0 0 0 0 0 1 11
12 22309.7 2.02 0 0 0 0 0 0 0 0 0 0 0 12
13 21786.3 2.03 1 0 0 0 0 0 0 0 0 0 0 13
14 21927.6 2.05 0 1 0 0 0 0 0 0 0 0 0 14
15 20957.9 2.08 0 0 1 0 0 0 0 0 0 0 0 15
16 19726.0 2.07 0 0 0 1 0 0 0 0 0 0 0 16
17 21315.7 2.06 0 0 0 0 1 0 0 0 0 0 0 17
18 24771.5 2.05 0 0 0 0 0 1 0 0 0 0 0 18
19 22592.4 2.05 0 0 0 0 0 0 1 0 0 0 0 19
20 21942.1 2.05 0 0 0 0 0 0 0 1 0 0 0 20
21 23973.7 2.05 0 0 0 0 0 0 0 0 1 0 0 21
22 20815.7 2.05 0 0 0 0 0 0 0 0 0 1 0 22
23 19931.4 2.06 0 0 0 0 0 0 0 0 0 0 1 23
24 24436.8 2.06 0 0 0 0 0 0 0 0 0 0 0 24
25 22838.7 2.07 1 0 0 0 0 0 0 0 0 0 0 25
26 24465.3 2.07 0 1 0 0 0 0 0 0 0 0 0 26
27 23007.3 2.30 0 0 1 0 0 0 0 0 0 0 0 27
28 22720.8 2.31 0 0 0 1 0 0 0 0 0 0 0 28
29 23045.7 2.31 0 0 0 0 1 0 0 0 0 0 0 29
30 27198.5 2.53 0 0 0 0 0 1 0 0 0 0 0 30
31 22401.9 2.58 0 0 0 0 0 0 1 0 0 0 0 31
32 25122.7 2.59 0 0 0 0 0 0 0 1 0 0 0 32
33 26100.5 2.73 0 0 0 0 0 0 0 0 1 0 0 33
34 22904.9 2.82 0 0 0 0 0 0 0 0 0 1 0 34
35 22040.4 3.00 0 0 0 0 0 0 0 0 0 0 1 35
36 25981.5 3.04 0 0 0 0 0 0 0 0 0 0 0 36
37 26157.1 3.23 1 0 0 0 0 0 0 0 0 0 0 37
38 25975.4 3.32 0 1 0 0 0 0 0 0 0 0 0 38
39 22589.8 3.49 0 0 1 0 0 0 0 0 0 0 0 39
40 25370.4 3.57 0 0 0 1 0 0 0 0 0 0 0 40
41 25091.1 3.56 0 0 0 0 1 0 0 0 0 0 0 41
42 28760.9 3.72 0 0 0 0 0 1 0 0 0 0 0 42
43 24325.9 3.82 0 0 0 0 0 0 1 0 0 0 0 43
44 25821.7 3.82 0 0 0 0 0 0 0 1 0 0 0 44
45 27645.7 3.98 0 0 0 0 0 0 0 0 1 0 0 45
46 26296.9 4.06 0 0 0 0 0 0 0 0 0 1 0 46
47 24141.5 4.08 0 0 0 0 0 0 0 0 0 0 1 47
48 27268.1 4.19 0 0 0 0 0 0 0 0 0 0 0 48
49 29060.3 4.16 1 0 0 0 0 0 0 0 0 0 0 49
50 28226.4 4.17 0 1 0 0 0 0 0 0 0 0 0 50
51 23268.5 4.21 0 0 1 0 0 0 0 0 0 0 0 51
52 26938.2 4.21 0 0 0 1 0 0 0 0 0 0 0 52
53 27217.5 4.17 0 0 0 0 1 0 0 0 0 0 0 53
54 27540.5 4.19 0 0 0 0 0 1 0 0 0 0 0 54
55 29167.6 4.25 0 0 0 0 0 0 1 0 0 0 0 55
56 26671.5 4.25 0 0 0 0 0 0 0 1 0 0 0 56
57 30184.0 4.20 0 0 0 0 0 0 0 0 1 0 0 57
58 28422.3 4.33 0 0 0 0 0 0 0 0 0 1 0 58
59 23774.3 4.41 0 0 0 0 0 0 0 0 0 0 1 59
60 29601.0 4.56 0 0 0 0 0 0 0 0 0 0 0 60
61 28523.6 5.18 1 0 0 0 0 0 0 0 0 0 0 61
62 23622.0 3.42 0 1 0 0 0 0 0 0 0 0 0 62
63 21320.3 2.71 0 0 1 0 0 0 0 0 0 0 0 63
64 20423.6 2.29 0 0 0 1 0 0 0 0 0 0 0 64
65 21174.9 2.00 0 0 0 0 1 0 0 0 0 0 0 65
66 23050.2 1.64 0 0 0 0 0 1 0 0 0 0 0 66
67 21202.9 1.30 0 0 0 0 0 0 1 0 0 0 0 67
68 20476.4 1.08 0 0 0 0 0 0 0 1 0 0 0 68
69 23173.3 1.00 0 0 0 0 0 0 0 0 1 0 0 69
70 22468.0 1.00 0 0 0 0 0 0 0 0 0 1 0 70
71 19842.7 1.00 0 0 0 0 0 0 0 0 0 0 1 71
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
18289.78 2000.13 -805.03 -1293.75 -3468.54 -2766.42
M5 M6 M7 M8 M9 M10
-2025.77 716.71 -1595.71 -1697.86 636.58 -1675.08
M11 t
-4123.36 35.59
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2421.02 -929.36 23.34 1005.33 2624.84
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 18289.78 781.62 23.400 < 2e-16 ***
X 2000.13 170.62 11.723 < 2e-16 ***
M1 -805.03 813.59 -0.989 0.326610
M2 -1293.75 814.30 -1.589 0.117638
M3 -3468.54 814.62 -4.258 7.80e-05 ***
M4 -2766.42 815.41 -3.393 0.001264 **
M5 -2025.77 816.55 -2.481 0.016081 *
M6 716.71 816.78 0.877 0.383908
M7 -1595.71 817.59 -1.952 0.055889 .
M8 -1697.86 818.80 -2.074 0.042647 *
M9 636.58 818.84 0.777 0.440126
M10 -1675.08 818.52 -2.046 0.045332 *
M11 -4123.36 818.31 -5.039 5.04e-06 ***
t 35.59 8.63 4.124 0.000122 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1342 on 57 degrees of freedom
Multiple R-squared: 0.8506, Adjusted R-squared: 0.8166
F-statistic: 24.97 on 13 and 57 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.1413416955 0.282683391 0.8586583
[2,] 0.0555255594 0.111051119 0.9444744
[3,] 0.0195993872 0.039198774 0.9804006
[4,] 0.0085328018 0.017065604 0.9914672
[5,] 0.0071871236 0.014374247 0.9928129
[6,] 0.0090810165 0.018162033 0.9909190
[7,] 0.0054881347 0.010976269 0.9945119
[8,] 0.0020999608 0.004199922 0.9979000
[9,] 0.0027441867 0.005488373 0.9972558
[10,] 0.0055861006 0.011172201 0.9944139
[11,] 0.0042540694 0.008508139 0.9957459
[12,] 0.0019771601 0.003954320 0.9980228
[13,] 0.0013045794 0.002609159 0.9986954
[14,] 0.0008217035 0.001643407 0.9991783
[15,] 0.0049031093 0.009806219 0.9950969
[16,] 0.0126529024 0.025305805 0.9873471
[17,] 0.0069465302 0.013893060 0.9930535
[18,] 0.0067043933 0.013408787 0.9932956
[19,] 0.0041096452 0.008219290 0.9958904
[20,] 0.0021671976 0.004334395 0.9978328
[21,] 0.0011579154 0.002315831 0.9988421
[22,] 0.0005939170 0.001187834 0.9994061
[23,] 0.0018619850 0.003723970 0.9981380
[24,] 0.0023199036 0.004639807 0.9976801
[25,] 0.0011324670 0.002264934 0.9988675
[26,] 0.0011165868 0.002233174 0.9988834
[27,] 0.0028898747 0.005779749 0.9971101
[28,] 0.0014150950 0.002830190 0.9985849
[29,] 0.0013329371 0.002665874 0.9986671
[30,] 0.0044847237 0.008969447 0.9955153
[31,] 0.0041420562 0.008284112 0.9958579
[32,] 0.0303072947 0.060614589 0.9696927
[33,] 0.0332197331 0.066439466 0.9667803
[34,] 0.0218033341 0.043606668 0.9781967
[35,] 0.2693961954 0.538792391 0.7306038
[36,] 0.1840410402 0.368082080 0.8159590
[37,] 0.1065589651 0.213117930 0.8934410
[38,] 0.8086219097 0.382756181 0.1913781
> postscript(file="/var/www/html/rcomp/tmp/1hup41258476422.ps",horizontal=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/2ultp1258476422.ps",horizontal=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/3sx641258476422.ps",horizontal=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/4avxo1258476422.ps",horizontal=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/5w6ix1258476422.ps",horizontal=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 = 71
Frequency = 1
1 2 3 4 5 6
-1016.00961 -2412.57707 -515.68349 -1542.38992 -1027.12879 -553.49479
7 8 9 10 11 12
-566.76434 -1663.30233 -515.93189 -1126.85858 -1432.76993 -447.42140
13 14 15 16 17 18
-221.37953 333.04769 1442.53860 -507.06782 326.39197 1024.12730
19 20 21 22 23 24
1121.85776 538.11977 199.69020 -682.23648 826.15217 1172.60070
25 26 27 28 29 30
323.94256 2403.67244 2624.83677 1580.62768 1129.28615 2063.99090
31 32 33 34 35 36
-555.78530 2211.57539 539.32721 -560.21144 627.95461 330.09782
37 38 39 40 41 42
895.11575 986.53367 -599.89403 1282.98757 247.44737 819.16010
43 44 45 46 47 48
-1539.02275 23.33927 -842.71157 -75.44889 141.83843 -1110.52767
49 50 51 52 53 54
1511.11952 1110.34807 -1788.36235 1143.62990 726.69368 -1768.37498
55 56 57 58 59 60
2015.54750 -413.99049 828.48659 1082.84262 -1312.47803 55.25055
61 62 63 64 65 66
-1492.78868 -2421.02480 -1163.43550 -1957.78741 -1402.69039 -1585.40853
67 68 69 70 71
-475.83287 -695.74161 -208.86054 1361.91278 1149.30276
> postscript(file="/var/www/html/rcomp/tmp/6ju8h1258476422.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 71
Frequency = 1
lag(myerror, k = 1) myerror
0 -1016.00961 NA
1 -2412.57707 -1016.00961
2 -515.68349 -2412.57707
3 -1542.38992 -515.68349
4 -1027.12879 -1542.38992
5 -553.49479 -1027.12879
6 -566.76434 -553.49479
7 -1663.30233 -566.76434
8 -515.93189 -1663.30233
9 -1126.85858 -515.93189
10 -1432.76993 -1126.85858
11 -447.42140 -1432.76993
12 -221.37953 -447.42140
13 333.04769 -221.37953
14 1442.53860 333.04769
15 -507.06782 1442.53860
16 326.39197 -507.06782
17 1024.12730 326.39197
18 1121.85776 1024.12730
19 538.11977 1121.85776
20 199.69020 538.11977
21 -682.23648 199.69020
22 826.15217 -682.23648
23 1172.60070 826.15217
24 323.94256 1172.60070
25 2403.67244 323.94256
26 2624.83677 2403.67244
27 1580.62768 2624.83677
28 1129.28615 1580.62768
29 2063.99090 1129.28615
30 -555.78530 2063.99090
31 2211.57539 -555.78530
32 539.32721 2211.57539
33 -560.21144 539.32721
34 627.95461 -560.21144
35 330.09782 627.95461
36 895.11575 330.09782
37 986.53367 895.11575
38 -599.89403 986.53367
39 1282.98757 -599.89403
40 247.44737 1282.98757
41 819.16010 247.44737
42 -1539.02275 819.16010
43 23.33927 -1539.02275
44 -842.71157 23.33927
45 -75.44889 -842.71157
46 141.83843 -75.44889
47 -1110.52767 141.83843
48 1511.11952 -1110.52767
49 1110.34807 1511.11952
50 -1788.36235 1110.34807
51 1143.62990 -1788.36235
52 726.69368 1143.62990
53 -1768.37498 726.69368
54 2015.54750 -1768.37498
55 -413.99049 2015.54750
56 828.48659 -413.99049
57 1082.84262 828.48659
58 -1312.47803 1082.84262
59 55.25055 -1312.47803
60 -1492.78868 55.25055
61 -2421.02480 -1492.78868
62 -1163.43550 -2421.02480
63 -1957.78741 -1163.43550
64 -1402.69039 -1957.78741
65 -1585.40853 -1402.69039
66 -475.83287 -1585.40853
67 -695.74161 -475.83287
68 -208.86054 -695.74161
69 1361.91278 -208.86054
70 1149.30276 1361.91278
71 NA 1149.30276
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -2412.57707 -1016.00961
[2,] -515.68349 -2412.57707
[3,] -1542.38992 -515.68349
[4,] -1027.12879 -1542.38992
[5,] -553.49479 -1027.12879
[6,] -566.76434 -553.49479
[7,] -1663.30233 -566.76434
[8,] -515.93189 -1663.30233
[9,] -1126.85858 -515.93189
[10,] -1432.76993 -1126.85858
[11,] -447.42140 -1432.76993
[12,] -221.37953 -447.42140
[13,] 333.04769 -221.37953
[14,] 1442.53860 333.04769
[15,] -507.06782 1442.53860
[16,] 326.39197 -507.06782
[17,] 1024.12730 326.39197
[18,] 1121.85776 1024.12730
[19,] 538.11977 1121.85776
[20,] 199.69020 538.11977
[21,] -682.23648 199.69020
[22,] 826.15217 -682.23648
[23,] 1172.60070 826.15217
[24,] 323.94256 1172.60070
[25,] 2403.67244 323.94256
[26,] 2624.83677 2403.67244
[27,] 1580.62768 2624.83677
[28,] 1129.28615 1580.62768
[29,] 2063.99090 1129.28615
[30,] -555.78530 2063.99090
[31,] 2211.57539 -555.78530
[32,] 539.32721 2211.57539
[33,] -560.21144 539.32721
[34,] 627.95461 -560.21144
[35,] 330.09782 627.95461
[36,] 895.11575 330.09782
[37,] 986.53367 895.11575
[38,] -599.89403 986.53367
[39,] 1282.98757 -599.89403
[40,] 247.44737 1282.98757
[41,] 819.16010 247.44737
[42,] -1539.02275 819.16010
[43,] 23.33927 -1539.02275
[44,] -842.71157 23.33927
[45,] -75.44889 -842.71157
[46,] 141.83843 -75.44889
[47,] -1110.52767 141.83843
[48,] 1511.11952 -1110.52767
[49,] 1110.34807 1511.11952
[50,] -1788.36235 1110.34807
[51,] 1143.62990 -1788.36235
[52,] 726.69368 1143.62990
[53,] -1768.37498 726.69368
[54,] 2015.54750 -1768.37498
[55,] -413.99049 2015.54750
[56,] 828.48659 -413.99049
[57,] 1082.84262 828.48659
[58,] -1312.47803 1082.84262
[59,] 55.25055 -1312.47803
[60,] -1492.78868 55.25055
[61,] -2421.02480 -1492.78868
[62,] -1163.43550 -2421.02480
[63,] -1957.78741 -1163.43550
[64,] -1402.69039 -1957.78741
[65,] -1585.40853 -1402.69039
[66,] -475.83287 -1585.40853
[67,] -695.74161 -475.83287
[68,] -208.86054 -695.74161
[69,] 1361.91278 -208.86054
[70,] 1149.30276 1361.91278
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -2412.57707 -1016.00961
2 -515.68349 -2412.57707
3 -1542.38992 -515.68349
4 -1027.12879 -1542.38992
5 -553.49479 -1027.12879
6 -566.76434 -553.49479
7 -1663.30233 -566.76434
8 -515.93189 -1663.30233
9 -1126.85858 -515.93189
10 -1432.76993 -1126.85858
11 -447.42140 -1432.76993
12 -221.37953 -447.42140
13 333.04769 -221.37953
14 1442.53860 333.04769
15 -507.06782 1442.53860
16 326.39197 -507.06782
17 1024.12730 326.39197
18 1121.85776 1024.12730
19 538.11977 1121.85776
20 199.69020 538.11977
21 -682.23648 199.69020
22 826.15217 -682.23648
23 1172.60070 826.15217
24 323.94256 1172.60070
25 2403.67244 323.94256
26 2624.83677 2403.67244
27 1580.62768 2624.83677
28 1129.28615 1580.62768
29 2063.99090 1129.28615
30 -555.78530 2063.99090
31 2211.57539 -555.78530
32 539.32721 2211.57539
33 -560.21144 539.32721
34 627.95461 -560.21144
35 330.09782 627.95461
36 895.11575 330.09782
37 986.53367 895.11575
38 -599.89403 986.53367
39 1282.98757 -599.89403
40 247.44737 1282.98757
41 819.16010 247.44737
42 -1539.02275 819.16010
43 23.33927 -1539.02275
44 -842.71157 23.33927
45 -75.44889 -842.71157
46 141.83843 -75.44889
47 -1110.52767 141.83843
48 1511.11952 -1110.52767
49 1110.34807 1511.11952
50 -1788.36235 1110.34807
51 1143.62990 -1788.36235
52 726.69368 1143.62990
53 -1768.37498 726.69368
54 2015.54750 -1768.37498
55 -413.99049 2015.54750
56 828.48659 -413.99049
57 1082.84262 828.48659
58 -1312.47803 1082.84262
59 55.25055 -1312.47803
60 -1492.78868 55.25055
61 -2421.02480 -1492.78868
62 -1163.43550 -2421.02480
63 -1957.78741 -1163.43550
64 -1402.69039 -1957.78741
65 -1585.40853 -1402.69039
66 -475.83287 -1585.40853
67 -695.74161 -475.83287
68 -208.86054 -695.74161
69 1361.91278 -208.86054
70 1149.30276 1361.91278
> 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/7mhra1258476422.ps",horizontal=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/8dsdc1258476422.ps",horizontal=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/940ye1258476422.ps",horizontal=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/10x6kl1258476422.ps",horizontal=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/11y4pq1258476422.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/12k4j01258476422.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/13ubdt1258476423.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/1421yg1258476423.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/158ciz1258476423.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/16y4lg1258476423.tab")
+ }
>
> system("convert tmp/1hup41258476422.ps tmp/1hup41258476422.png")
> system("convert tmp/2ultp1258476422.ps tmp/2ultp1258476422.png")
> system("convert tmp/3sx641258476422.ps tmp/3sx641258476422.png")
> system("convert tmp/4avxo1258476422.ps tmp/4avxo1258476422.png")
> system("convert tmp/5w6ix1258476422.ps tmp/5w6ix1258476422.png")
> system("convert tmp/6ju8h1258476422.ps tmp/6ju8h1258476422.png")
> system("convert tmp/7mhra1258476422.ps tmp/7mhra1258476422.png")
> system("convert tmp/8dsdc1258476422.ps tmp/8dsdc1258476422.png")
> system("convert tmp/940ye1258476422.ps tmp/940ye1258476422.png")
> system("convert tmp/10x6kl1258476422.ps tmp/10x6kl1258476422.png")
>
>
> proc.time()
user system elapsed
2.590 1.598 3.506