R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-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(2050
+ ,2650
+ ,13
+ ,7
+ ,1
+ ,1
+ ,0
+ ,1639
+ ,2150
+ ,2664
+ ,6
+ ,5
+ ,1
+ ,1
+ ,0
+ ,1193
+ ,2150
+ ,2921
+ ,3
+ ,6
+ ,1
+ ,1
+ ,0
+ ,1635
+ ,1999
+ ,2580
+ ,4
+ ,4
+ ,1
+ ,1
+ ,0
+ ,1732
+ ,1900
+ ,2580
+ ,4
+ ,4
+ ,1
+ ,0
+ ,0
+ ,1534
+ ,1800
+ ,2774
+ ,2
+ ,4
+ ,1
+ ,0
+ ,0
+ ,1765
+ ,1560
+ ,1920
+ ,1
+ ,5
+ ,1
+ ,1
+ ,0
+ ,1161
+ ,1449
+ ,1710
+ ,1
+ ,3
+ ,1
+ ,1
+ ,0
+ ,1010
+ ,1375
+ ,1837
+ ,4
+ ,5
+ ,1
+ ,0
+ ,0
+ ,1191
+ ,1270
+ ,1880
+ ,8
+ ,6
+ ,1
+ ,0
+ ,0
+ ,930
+ ,1250
+ ,2150
+ ,15
+ ,3
+ ,1
+ ,0
+ ,0
+ ,984
+ ,1235
+ ,1894
+ ,14
+ ,5
+ ,1
+ ,1
+ ,0
+ ,1112
+ ,1170
+ ,1928
+ ,18
+ ,8
+ ,1
+ ,1
+ ,0
+ ,600
+ ,1155
+ ,1767
+ ,16
+ ,4
+ ,1
+ ,0
+ ,0
+ ,794
+ ,1110
+ ,1630
+ ,15
+ ,3
+ ,1
+ ,0
+ ,1
+ ,867
+ ,1139
+ ,1680
+ ,17
+ ,4
+ ,1
+ ,0
+ ,1
+ ,750
+ ,995
+ ,1500
+ ,15
+ ,4
+ ,1
+ ,0
+ ,0
+ ,743
+ ,900
+ ,1400
+ ,16
+ ,2
+ ,1
+ ,0
+ ,1
+ ,731
+ ,960
+ ,1573
+ ,17
+ ,6
+ ,1
+ ,0
+ ,0
+ ,768
+ ,1695
+ ,2931
+ ,28
+ ,3
+ ,1
+ ,0
+ ,1
+ ,1142
+ ,1553
+ ,2200
+ ,28
+ ,4
+ ,1
+ ,0
+ ,0
+ ,1035
+ ,1020
+ ,1478
+ ,53
+ ,3
+ ,1
+ ,0
+ ,1
+ ,626
+ ,1020
+ ,1713
+ ,30
+ ,4
+ ,1
+ ,0
+ ,1
+ ,600
+ ,850
+ ,1190
+ ,41
+ ,1
+ ,1
+ ,0
+ ,0
+ ,600
+ ,720
+ ,1121
+ ,46
+ ,4
+ ,1
+ ,0
+ ,0
+ ,398
+ ,749
+ ,1733
+ ,43
+ ,6
+ ,1
+ ,0
+ ,0
+ ,656
+ ,2150
+ ,2848
+ ,4
+ ,6
+ ,1
+ ,1
+ ,0
+ ,1487
+ ,1350
+ ,2253
+ ,23
+ ,4
+ ,1
+ ,1
+ ,0
+ ,939
+ ,1299
+ ,2743
+ ,25
+ ,5
+ ,1
+ ,1
+ ,1
+ ,1232
+ ,1250
+ ,2180
+ ,17
+ ,4
+ ,1
+ ,0
+ ,1
+ ,1141
+ ,1239
+ ,1706
+ ,14
+ ,4
+ ,1
+ ,0
+ ,0
+ ,810
+ ,1125
+ ,1710
+ ,16
+ ,4
+ ,1
+ ,1
+ ,0
+ ,800
+ ,1080
+ ,2200
+ ,26
+ ,4
+ ,1
+ ,0
+ ,0
+ ,1076
+ ,1050
+ ,1680
+ ,13
+ ,4
+ ,1
+ ,0
+ ,0
+ ,875
+ ,1049
+ ,1900
+ ,34
+ ,3
+ ,1
+ ,0
+ ,0
+ ,690
+ ,934
+ ,1543
+ ,20
+ ,3
+ ,1
+ ,0
+ ,0
+ ,820
+ ,875
+ ,1173
+ ,6
+ ,4
+ ,1
+ ,0
+ ,0
+ ,456
+ ,805
+ ,1258
+ ,7
+ ,4
+ ,1
+ ,0
+ ,1
+ ,821
+ ,759
+ ,997
+ ,4
+ ,4
+ ,1
+ ,0
+ ,0
+ ,461
+ ,729
+ ,1007
+ ,19
+ ,6
+ ,1
+ ,0
+ ,0
+ ,513
+ ,710
+ ,1083
+ ,22
+ ,4
+ ,1
+ ,0
+ ,0
+ ,504
+ ,975
+ ,1500
+ ,7
+ ,3
+ ,0
+ ,1
+ ,1
+ ,700
+ ,939
+ ,1428
+ ,40
+ ,2
+ ,0
+ ,0
+ ,0
+ ,701
+ ,2100
+ ,2116
+ ,25
+ ,3
+ ,0
+ ,1
+ ,0
+ ,1209
+ ,580
+ ,1051
+ ,15
+ ,2
+ ,0
+ ,0
+ ,0
+ ,426
+ ,1844
+ ,2250
+ ,40
+ ,6
+ ,0
+ ,1
+ ,0
+ ,915
+ ,699
+ ,1400
+ ,45
+ ,1
+ ,0
+ ,1
+ ,1
+ ,481
+ ,1160
+ ,1720
+ ,5
+ ,4
+ ,0
+ ,0
+ ,0
+ ,867
+ ,1109
+ ,1740
+ ,4
+ ,3
+ ,0
+ ,0
+ ,0
+ ,816
+ ,1129
+ ,1700
+ ,6
+ ,4
+ ,0
+ ,0
+ ,0
+ ,725
+ ,1050
+ ,1620
+ ,6
+ ,4
+ ,0
+ ,0
+ ,0
+ ,800
+ ,1045
+ ,1630
+ ,6
+ ,4
+ ,0
+ ,0
+ ,0
+ ,750
+ ,1050
+ ,1920
+ ,8
+ ,4
+ ,0
+ ,0
+ ,0
+ ,944
+ ,1020
+ ,1606
+ ,5
+ ,4
+ ,0
+ ,0
+ ,0
+ ,811
+ ,1000
+ ,1535
+ ,7
+ ,5
+ ,0
+ ,0
+ ,1
+ ,668
+ ,1030
+ ,1540
+ ,6
+ ,2
+ ,0
+ ,0
+ ,1
+ ,826
+ ,975
+ ,1739
+ ,13
+ ,3
+ ,0
+ ,0
+ ,0
+ ,880
+ ,940
+ ,1305
+ ,5
+ ,3
+ ,0
+ ,0
+ ,0
+ ,647
+ ,920
+ ,1415
+ ,7
+ ,4
+ ,0
+ ,0
+ ,0
+ ,866
+ ,945
+ ,1580
+ ,9
+ ,3
+ ,0
+ ,0
+ ,0
+ ,810
+ ,874
+ ,1236
+ ,3
+ ,4
+ ,0
+ ,0
+ ,0
+ ,707
+ ,872
+ ,1229
+ ,6
+ ,3
+ ,0
+ ,0
+ ,0
+ ,721
+ ,870
+ ,1273
+ ,4
+ ,4
+ ,0
+ ,0
+ ,0
+ ,638
+ ,869
+ ,1165
+ ,7
+ ,4
+ ,0
+ ,0
+ ,0
+ ,694
+ ,766
+ ,1200
+ ,7
+ ,4
+ ,0
+ ,0
+ ,1
+ ,634
+ ,739
+ ,970
+ ,4
+ ,4
+ ,0
+ ,0
+ ,1
+ ,541)
+ ,dim=c(8
+ ,66)
+ ,dimnames=list(c('PRICE'
+ ,'SQFT'
+ ,'AGE'
+ ,'FEATS'
+ ,'NE'
+ ,'CUST'
+ ,'COR'
+ ,'TAX')
+ ,1:66))
> y <- array(NA,dim=c(8,66),dimnames=list(c('PRICE','SQFT','AGE','FEATS','NE','CUST','COR','TAX'),1:66))
> 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 = '1'
> par3 <- 'No 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, 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
PRICE SQFT AGE FEATS NE CUST COR TAX
1 2050 2650 13 7 1 1 0 1639
2 2150 2664 6 5 1 1 0 1193
3 2150 2921 3 6 1 1 0 1635
4 1999 2580 4 4 1 1 0 1732
5 1900 2580 4 4 1 0 0 1534
6 1800 2774 2 4 1 0 0 1765
7 1560 1920 1 5 1 1 0 1161
8 1449 1710 1 3 1 1 0 1010
9 1375 1837 4 5 1 0 0 1191
10 1270 1880 8 6 1 0 0 930
11 1250 2150 15 3 1 0 0 984
12 1235 1894 14 5 1 1 0 1112
13 1170 1928 18 8 1 1 0 600
14 1155 1767 16 4 1 0 0 794
15 1110 1630 15 3 1 0 1 867
16 1139 1680 17 4 1 0 1 750
17 995 1500 15 4 1 0 0 743
18 900 1400 16 2 1 0 1 731
19 960 1573 17 6 1 0 0 768
20 1695 2931 28 3 1 0 1 1142
21 1553 2200 28 4 1 0 0 1035
22 1020 1478 53 3 1 0 1 626
23 1020 1713 30 4 1 0 1 600
24 850 1190 41 1 1 0 0 600
25 720 1121 46 4 1 0 0 398
26 749 1733 43 6 1 0 0 656
27 2150 2848 4 6 1 1 0 1487
28 1350 2253 23 4 1 1 0 939
29 1299 2743 25 5 1 1 1 1232
30 1250 2180 17 4 1 0 1 1141
31 1239 1706 14 4 1 0 0 810
32 1125 1710 16 4 1 1 0 800
33 1080 2200 26 4 1 0 0 1076
34 1050 1680 13 4 1 0 0 875
35 1049 1900 34 3 1 0 0 690
36 934 1543 20 3 1 0 0 820
37 875 1173 6 4 1 0 0 456
38 805 1258 7 4 1 0 1 821
39 759 997 4 4 1 0 0 461
40 729 1007 19 6 1 0 0 513
41 710 1083 22 4 1 0 0 504
42 975 1500 7 3 0 1 1 700
43 939 1428 40 2 0 0 0 701
44 2100 2116 25 3 0 1 0 1209
45 580 1051 15 2 0 0 0 426
46 1844 2250 40 6 0 1 0 915
47 699 1400 45 1 0 1 1 481
48 1160 1720 5 4 0 0 0 867
49 1109 1740 4 3 0 0 0 816
50 1129 1700 6 4 0 0 0 725
51 1050 1620 6 4 0 0 0 800
52 1045 1630 6 4 0 0 0 750
53 1050 1920 8 4 0 0 0 944
54 1020 1606 5 4 0 0 0 811
55 1000 1535 7 5 0 0 1 668
56 1030 1540 6 2 0 0 1 826
57 975 1739 13 3 0 0 0 880
58 940 1305 5 3 0 0 0 647
59 920 1415 7 4 0 0 0 866
60 945 1580 9 3 0 0 0 810
61 874 1236 3 4 0 0 0 707
62 872 1229 6 3 0 0 0 721
63 870 1273 4 4 0 0 0 638
64 869 1165 7 4 0 0 0 694
65 766 1200 7 4 0 0 1 634
66 739 970 4 4 0 0 1 541
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) SQFT AGE FEATS NE CUST
92.7448 0.3522 -0.5651 4.3896 -17.3853 174.9411
COR TAX
-73.5823 0.4989
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-466.28 -82.29 6.75 78.70 484.84
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 92.74480 101.60704 0.913 0.365137
SQFT 0.35222 0.09575 3.679 0.000515 ***
AGE -0.56508 2.00253 -0.282 0.778807
FEATS 4.38961 18.55499 0.237 0.813822
NE -17.38534 47.27462 -0.368 0.714397
CUST 174.94108 53.72371 3.256 0.001887 **
COR -73.58234 49.13007 -1.498 0.139633
TAX 0.49887 0.15849 3.148 0.002598 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 158.9 on 58 degrees of freedom
Multiple R-squared: 0.8623, Adjusted R-squared: 0.8456
F-statistic: 51.86 on 7 and 58 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.57974342 8.405132e-01 4.202566e-01
[2,] 0.58028806 8.394239e-01 4.197119e-01
[3,] 0.55655986 8.868803e-01 4.434401e-01
[4,] 0.50573667 9.885267e-01 4.942633e-01
[5,] 0.38593969 7.718794e-01 6.140603e-01
[6,] 0.29615377 5.923075e-01 7.038462e-01
[7,] 0.23416121 4.683224e-01 7.658388e-01
[8,] 0.16583580 3.316716e-01 8.341642e-01
[9,] 0.11001897 2.200379e-01 8.899810e-01
[10,] 0.08390088 1.678018e-01 9.160991e-01
[11,] 0.23191535 4.638307e-01 7.680847e-01
[12,] 0.28398606 5.679721e-01 7.160139e-01
[13,] 0.27820611 5.564122e-01 7.217939e-01
[14,] 0.22499602 4.499920e-01 7.750040e-01
[15,] 0.17446007 3.489201e-01 8.255399e-01
[16,] 0.30146848 6.029370e-01 6.985315e-01
[17,] 0.27707770 5.541554e-01 7.229223e-01
[18,] 0.33135544 6.627109e-01 6.686446e-01
[19,] 0.84977285 3.004543e-01 1.502272e-01
[20,] 0.81323218 3.735356e-01 1.867678e-01
[21,] 0.81582271 3.683546e-01 1.841773e-01
[22,] 0.84999152 3.000170e-01 1.500085e-01
[23,] 0.94475789 1.104842e-01 5.524211e-02
[24,] 0.93117508 1.376498e-01 6.882492e-02
[25,] 0.90098142 1.980372e-01 9.901858e-02
[26,] 0.88285411 2.342918e-01 1.171459e-01
[27,] 0.88176121 2.364776e-01 1.182388e-01
[28,] 0.87001329 2.599734e-01 1.299867e-01
[29,] 0.85810293 2.837941e-01 1.418971e-01
[30,] 0.81411930 3.717614e-01 1.858807e-01
[31,] 0.75075537 4.984893e-01 2.492446e-01
[32,] 0.78859373 4.228125e-01 2.114063e-01
[33,] 0.97279824 5.440351e-02 2.720176e-02
[34,] 0.99335527 1.328946e-02 6.644731e-03
[35,] 0.99631950 7.361007e-03 3.680503e-03
[36,] 0.99998571 2.857277e-05 1.428639e-05
[37,] 0.99995799 8.402403e-05 4.201201e-05
[38,] 0.99995424 9.151967e-05 4.575984e-05
[39,] 0.99983514 3.297187e-04 1.648594e-04
[40,] 0.99964560 7.087969e-04 3.543985e-04
[41,] 0.99896882 2.062367e-03 1.031183e-03
[42,] 0.99673837 6.523252e-03 3.261626e-03
[43,] 0.99420576 1.158848e-02 5.794242e-03
[44,] 0.97997950 4.004101e-02 2.002050e-02
[45,] 0.98180810 3.638380e-02 1.819190e-02
> postscript(file="/var/wessaorg/rcomp/tmp/1ww3f1355130893.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/wessaorg/rcomp/tmp/2le431355130893.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/wessaorg/rcomp/tmp/30gsp1355130893.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/wessaorg/rcomp/tmp/44yyx1355130893.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/wessaorg/rcomp/tmp/5o8fk1355130893.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 = 66
Frequency = 1
1 2 3 4 5 6
25.291615 347.680244 30.574711 -39.364930 135.352420 -149.347086
7 8 9 10 11 12
32.869134 79.943580 38.773480 46.703895 -78.209653 -251.182495
13 14 15 16 17 18
-83.644945 52.650763 96.894031 163.391466 11.570355 40.705269
19 20 21 22 23 24
-69.262386 93.814747 184.693512 202.132007 114.944830 74.957406
25 26 27 28 29 30
59.688812 -266.051754 130.684501 -166.849024 -466.282044 -96.775902
31 32 33 34 35 36
149.023998 -135.207091 -309.890324 -63.809960 -33.750716 -95.773014
37 38 39 40 41 42
144.835727 -63.042978 87.201644 27.435230 -3.369076 -105.853354
43 44 45 46 47 48
7.403536 484.844409 -95.747921 323.622350 -207.126661 14.186275
49 50 51 52 53 54
-14.591194 61.635272 -26.602514 -10.181195 -202.975145 -57.724109
55 56 57 58 59 60
88.944711 50.965889 -175.080924 54.497911 -116.758093 -116.417628
61 62 63 64 65 66
-22.650995 -23.084796 -4.695959 6.102145 -5.710951 92.998940
> postscript(file="/var/wessaorg/rcomp/tmp/6w7jm1355130893.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 = 66
Frequency = 1
lag(myerror, k = 1) myerror
0 25.291615 NA
1 347.680244 25.291615
2 30.574711 347.680244
3 -39.364930 30.574711
4 135.352420 -39.364930
5 -149.347086 135.352420
6 32.869134 -149.347086
7 79.943580 32.869134
8 38.773480 79.943580
9 46.703895 38.773480
10 -78.209653 46.703895
11 -251.182495 -78.209653
12 -83.644945 -251.182495
13 52.650763 -83.644945
14 96.894031 52.650763
15 163.391466 96.894031
16 11.570355 163.391466
17 40.705269 11.570355
18 -69.262386 40.705269
19 93.814747 -69.262386
20 184.693512 93.814747
21 202.132007 184.693512
22 114.944830 202.132007
23 74.957406 114.944830
24 59.688812 74.957406
25 -266.051754 59.688812
26 130.684501 -266.051754
27 -166.849024 130.684501
28 -466.282044 -166.849024
29 -96.775902 -466.282044
30 149.023998 -96.775902
31 -135.207091 149.023998
32 -309.890324 -135.207091
33 -63.809960 -309.890324
34 -33.750716 -63.809960
35 -95.773014 -33.750716
36 144.835727 -95.773014
37 -63.042978 144.835727
38 87.201644 -63.042978
39 27.435230 87.201644
40 -3.369076 27.435230
41 -105.853354 -3.369076
42 7.403536 -105.853354
43 484.844409 7.403536
44 -95.747921 484.844409
45 323.622350 -95.747921
46 -207.126661 323.622350
47 14.186275 -207.126661
48 -14.591194 14.186275
49 61.635272 -14.591194
50 -26.602514 61.635272
51 -10.181195 -26.602514
52 -202.975145 -10.181195
53 -57.724109 -202.975145
54 88.944711 -57.724109
55 50.965889 88.944711
56 -175.080924 50.965889
57 54.497911 -175.080924
58 -116.758093 54.497911
59 -116.417628 -116.758093
60 -22.650995 -116.417628
61 -23.084796 -22.650995
62 -4.695959 -23.084796
63 6.102145 -4.695959
64 -5.710951 6.102145
65 92.998940 -5.710951
66 NA 92.998940
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 347.680244 25.291615
[2,] 30.574711 347.680244
[3,] -39.364930 30.574711
[4,] 135.352420 -39.364930
[5,] -149.347086 135.352420
[6,] 32.869134 -149.347086
[7,] 79.943580 32.869134
[8,] 38.773480 79.943580
[9,] 46.703895 38.773480
[10,] -78.209653 46.703895
[11,] -251.182495 -78.209653
[12,] -83.644945 -251.182495
[13,] 52.650763 -83.644945
[14,] 96.894031 52.650763
[15,] 163.391466 96.894031
[16,] 11.570355 163.391466
[17,] 40.705269 11.570355
[18,] -69.262386 40.705269
[19,] 93.814747 -69.262386
[20,] 184.693512 93.814747
[21,] 202.132007 184.693512
[22,] 114.944830 202.132007
[23,] 74.957406 114.944830
[24,] 59.688812 74.957406
[25,] -266.051754 59.688812
[26,] 130.684501 -266.051754
[27,] -166.849024 130.684501
[28,] -466.282044 -166.849024
[29,] -96.775902 -466.282044
[30,] 149.023998 -96.775902
[31,] -135.207091 149.023998
[32,] -309.890324 -135.207091
[33,] -63.809960 -309.890324
[34,] -33.750716 -63.809960
[35,] -95.773014 -33.750716
[36,] 144.835727 -95.773014
[37,] -63.042978 144.835727
[38,] 87.201644 -63.042978
[39,] 27.435230 87.201644
[40,] -3.369076 27.435230
[41,] -105.853354 -3.369076
[42,] 7.403536 -105.853354
[43,] 484.844409 7.403536
[44,] -95.747921 484.844409
[45,] 323.622350 -95.747921
[46,] -207.126661 323.622350
[47,] 14.186275 -207.126661
[48,] -14.591194 14.186275
[49,] 61.635272 -14.591194
[50,] -26.602514 61.635272
[51,] -10.181195 -26.602514
[52,] -202.975145 -10.181195
[53,] -57.724109 -202.975145
[54,] 88.944711 -57.724109
[55,] 50.965889 88.944711
[56,] -175.080924 50.965889
[57,] 54.497911 -175.080924
[58,] -116.758093 54.497911
[59,] -116.417628 -116.758093
[60,] -22.650995 -116.417628
[61,] -23.084796 -22.650995
[62,] -4.695959 -23.084796
[63,] 6.102145 -4.695959
[64,] -5.710951 6.102145
[65,] 92.998940 -5.710951
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 347.680244 25.291615
2 30.574711 347.680244
3 -39.364930 30.574711
4 135.352420 -39.364930
5 -149.347086 135.352420
6 32.869134 -149.347086
7 79.943580 32.869134
8 38.773480 79.943580
9 46.703895 38.773480
10 -78.209653 46.703895
11 -251.182495 -78.209653
12 -83.644945 -251.182495
13 52.650763 -83.644945
14 96.894031 52.650763
15 163.391466 96.894031
16 11.570355 163.391466
17 40.705269 11.570355
18 -69.262386 40.705269
19 93.814747 -69.262386
20 184.693512 93.814747
21 202.132007 184.693512
22 114.944830 202.132007
23 74.957406 114.944830
24 59.688812 74.957406
25 -266.051754 59.688812
26 130.684501 -266.051754
27 -166.849024 130.684501
28 -466.282044 -166.849024
29 -96.775902 -466.282044
30 149.023998 -96.775902
31 -135.207091 149.023998
32 -309.890324 -135.207091
33 -63.809960 -309.890324
34 -33.750716 -63.809960
35 -95.773014 -33.750716
36 144.835727 -95.773014
37 -63.042978 144.835727
38 87.201644 -63.042978
39 27.435230 87.201644
40 -3.369076 27.435230
41 -105.853354 -3.369076
42 7.403536 -105.853354
43 484.844409 7.403536
44 -95.747921 484.844409
45 323.622350 -95.747921
46 -207.126661 323.622350
47 14.186275 -207.126661
48 -14.591194 14.186275
49 61.635272 -14.591194
50 -26.602514 61.635272
51 -10.181195 -26.602514
52 -202.975145 -10.181195
53 -57.724109 -202.975145
54 88.944711 -57.724109
55 50.965889 88.944711
56 -175.080924 50.965889
57 54.497911 -175.080924
58 -116.758093 54.497911
59 -116.417628 -116.758093
60 -22.650995 -116.417628
61 -23.084796 -22.650995
62 -4.695959 -23.084796
63 6.102145 -4.695959
64 -5.710951 6.102145
65 92.998940 -5.710951
> 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/wessaorg/rcomp/tmp/7g3261355130893.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/wessaorg/rcomp/tmp/82nfk1355130893.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/wessaorg/rcomp/tmp/9cvw31355130893.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/wessaorg/rcomp/tmp/10agss1355130893.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/1181lz1355130893.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/wessaorg/rcomp/tmp/12r3h81355130894.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/wessaorg/rcomp/tmp/13ek5k1355130894.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/wessaorg/rcomp/tmp/147tgd1355130894.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/wessaorg/rcomp/tmp/15fq9c1355130894.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/wessaorg/rcomp/tmp/16750q1355130894.tab")
+ }
>
> try(system("convert tmp/1ww3f1355130893.ps tmp/1ww3f1355130893.png",intern=TRUE))
character(0)
> try(system("convert tmp/2le431355130893.ps tmp/2le431355130893.png",intern=TRUE))
character(0)
> try(system("convert tmp/30gsp1355130893.ps tmp/30gsp1355130893.png",intern=TRUE))
character(0)
> try(system("convert tmp/44yyx1355130893.ps tmp/44yyx1355130893.png",intern=TRUE))
character(0)
> try(system("convert tmp/5o8fk1355130893.ps tmp/5o8fk1355130893.png",intern=TRUE))
character(0)
> try(system("convert tmp/6w7jm1355130893.ps tmp/6w7jm1355130893.png",intern=TRUE))
character(0)
> try(system("convert tmp/7g3261355130893.ps tmp/7g3261355130893.png",intern=TRUE))
character(0)
> try(system("convert tmp/82nfk1355130893.ps tmp/82nfk1355130893.png",intern=TRUE))
character(0)
> try(system("convert tmp/9cvw31355130893.ps tmp/9cvw31355130893.png",intern=TRUE))
character(0)
> try(system("convert tmp/10agss1355130893.ps tmp/10agss1355130893.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
5.967 0.877 6.838