R语言代码试题答案步骤.doc

R version 3.4.3 (2017-11-30) -- "Kite-Eating Tree"
Copyright (C) 2017 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R是自由软件，不带任何担保。
在某些条件下你可以将其自由散布。
用'license()'或'licence()'来看散布的详细条件。

R是个合作计划，有许多人为之做出了贡献.
用'contributors()'来看合作者的详细情况
用'citation()'会告诉你如何在出版物中正确地引用R或R程序包。

用'demo()'来看一些示范程序，用'help()'来阅读在线帮助文件，或
用'help.start()'通过HTML浏览器来看帮助文件。
用'q()'退出R.

[原来保存的工作空间已还原]

> h=read.csv("F://1.csv",header=true)
Error in read.table(file = file, header = header, sep = sep, quote = quote, :
找不到对象'true'
> h=read.csv("F://1.csv",header=TRUE)
> h
地区 x1 x2 x3 x4 x5 x6 x7 x8 x9 y
1 北京 7535 2639 1971 1658 3696 84742 87475 106.5 1.3 24046
2 天津 7344 1881 1854 1556 2254 61514 93173 107.5 3.6 20024
3 河北 4211 1542 1502 1047 1204 38658 36584 104.1 3.7 12531
4 山西 3856 1529 1439 906 1506 44236 33628 108.8 3.3 12212
5 内蒙古 5463 2730 1584 1354 1972 46557 63886 109.6 3.7 17717
6 辽宁 5809 2042 1433 1310 1844 41858 56649 107.7 3.6 16594
7 吉林 4635 2045 1594 1448 1643 38407 43415 111.0 3.7 14614
8 黑龙江 4687 1807 1337 1181 1217 36406 35711 104.8 4.2 12984
9 上海 9656 2111 1790 1017 3724 78673 85373 106.0 3.1 26253
10 江苏 6658 1916 1437 1058 3078 50639 68347 112.6 3.1 18825
11 浙江 7552 2110 1552 1228 2997 50197 63374 104.5 3.0 21545
12 安徽 5815 1541 1397 1143 1933 44601 28792 105.3 3.7 15012
13 福建 7317 1634 1754 773 2105 44525 52763 104.6 3.6 18593
14 江西 5072 1477 1174 671 1487 38512 28800 106.7 3.0 12776
15 山东 5201 2197 1572 1005 1656 41904 51768 106.9 3.3 15778
16 河南 4607 1886 1191 1085 1525 37338 31499 106.8 3.1 13733
17 湖北 5838 1783 1371 1030 1652 39846 38572 105.6 3.8 14496
18 湖南 5442 1625 1302 918 1738 38971 33480 105.7 4.2 14609
19 广东 8258 1521 2100 1048 2954 50278 54095 107.9 2.5 22396
20 广西 5553 1146 1377 884 1626 36386 27952 107.5 3.4 14244
21 海南 6556 865 1521 993 1320 39485 32377 107.0 2.0 14457
22 重庆 6870 2229 1177 1102 1471 44498 38914 107.8 3.3 16573
23 四川 6074 1651 1284 773 1587 42339 29608 105.9 4.0 15050
24 贵州 4993 1399 1014 655 1396 41156 19710 105.5 3.3 12586
25 云南 5468 1760 974 939 1434 37629 22195 108.9 4.0 13884
26 西藏 5518 1362 845 467 550 51705 22936 109.5 2.6 11184
27 陕西 5551 1789 1322 1212 2079 43073 38564 109.4 3.2 15333
28 甘肃 4602 1631 1288 1050 1388 37679 21978 108.6 2.7 12847
29 青海 4667 1512 1232 906 1097 46483 33181 110.6 3.4 12346
30 宁夏 4769 1876 1193 1063 1516 47436 36394 105.5 4.2 14067
31 新疆 5239 2031 1167 1028 1281 44576 33796 114.8 3.4 13892
> lm=lm(y~x1+x2+x3+x4+x5+x6+x7+x8+x9,data=h)
> lm

Call:
lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9,
data = h)

Coefficients:
(Intercept) x1 x2 x3 x4 x5 x6 x7 x8 x9
320.640948 1.316588 1.649859 2.178660 -0.005609 1.684283 0.010320 0.003655 -19.130576 50.515575

> summary(lm)

Call:
lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9,
data = h)

Residuals:
Min 1Q Median 3Q Max
-940.13 -195.24 3.42 239.00 476.06

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.206e+02 3.952e+03 0.081 0.936097
x1 1.317e+00 1.062e-01 12.400 3.97e-11 ***
x2 1.650e+00 3.008e-01 5.484 1.93e-05 ***
x3 2.179e+00 5.199e-01 4.190 0.000412 ***
x4 -5.609e-03 4.766e-01 -0.012 0.990720
x5 1.684e+00 2.142e-01 7.864 1.08e-07 ***
x6 1.032e-02 1.343e-02 0.769 0.450665
x7 3.655e-03 1.070e-02 0.342 0.736006
x8 -1.913e+01 3.197e+01 -0.598 0.555983
x9 5.052e+01 1.502e+02 0.336 0.739986
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 389.4 on 21 degrees of freedom
Multiple R-squared: 0.9923, Adjusted R-squared: 0.9889
F-statistic: 298.9 on 9 and 21 DF, p-value: < 2.2e-16

> pre=fitted.values(lm)
> res=residuals(lm)
> sd(res)
[1] 325.7967
> res=residuals(lm)

> dy=step(lm)
Start: AIC=377.73
y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9

Df Sum of Sq RSS AIC
- x4 1 21 3184326 375.73
- x9 1 17149 3201454 375.90
- x7 1 17700 3202005 375.90
- x8 1 54295 3238599 376.26
- x6 1 89586 3273891 376.59
<none> 3184305 377.73
- x3 1 2662593 5846898 394.57
- x2 1 4561056 7745361 403.29
- x5 1 9377500 12561805 418.28
- x1 1 23314547 26498852 441.42

Step: AIC=375.73
y ~ x1 + x2 + x3 + x5 + x6 + x7 + x8 + x9

Df Sum of Sq RSS AIC
- x9 1 17428 3201754 373.90
- x7 1 18563 3202889 373.91
- x8 1 54437 3238763 374.26
- x6 1 91813 3276139 374.61
<none> 3184326 375.73
- x3 1 2936130 6120456 393.99
- x2 1 5467941 8652267 404.72
- x5 1 9393345 12577671 416.32
- x1 1 25886086 29070412 442.29

Step: AIC=373.9
y ~ x1 + x2 + x3 + x5 + x6 + x7 + x8

Df Sum of Sq RSS AIC
- x7 1 34634 3236387 372.24
- x6 1 74800 3276554 372.62
- x8 1 82150 3283904 372.69
<none> 3201754 373.90
- x3 1 3055353 6257107 392.67
- x2 1 5725836 8927590 403.69
- x5 1 9382624 12584378 414.33
- x1 1 25868832 29070586 440.29

Step: AIC=372.24
y ~ x1 + x2 + x3 + x5 + x6 + x8

Df Sum of Sq RSS AIC
- x8 1 70813 3307201 370.91
- x6 1 152777 3389165 371.67
<none> 3236387 372.24
- x3 1 5501284 8737672 401.02
- x2 1 8895049 12131436 411.20
- x5 1 9458098 12694485 412.60
- x1 1 27733098 30969486 440.25

Step: AIC=370.91
y ~ x1 + x2 + x3 + x5 + x6

Df Sum of Sq RSS AIC
- x6 1 137540 3444741 370.17
<none> 3307201 370.91
- x3 1 5771063 9078264 400.21
- x2 1 8871193 12178394 409.32
- x5 1 9473521 12780722 410.81
- x1 1 28248162 31555363 438.83

Step: AIC=370.17
y ~ x1 + x2 + x3 + x5

Df Sum of Sq RSS AIC
<none> 3444741 370.17
- x3 1 5717883 9162624 398.50
- x2 1 10249815 13694556 410.95
- x5 1 10998313 14443054 412.60
- x1 1 33258637 36703378 441.52
> summary(dy)

Call:
lm(formula = y ~ x1 + x2 + x3 + x5, data = h)

Residuals:
Min 1Q Median 3Q Max
-943.18 -161.05 12.74 250.93 566.25

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1694.6269 562.9773 -3.010 0.00574 **
x1 1.3642 0.0861 15.844 7.11e-15 ***
x2 1.7679 0.2010 8.796 2.86e-09 ***
x3 2.2894 0.3485 6.569 5.76e-07 ***
x5 1.7424 0.1912 9.111 1.42e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 364 on 26 degrees of freedom
Multiple R-squared: 0.9916, Adjusted R-squared: 0.9903
F-statistic: 769.2 on 4 and 26 DF, p-value: < 2.2e-16

> newdata=data.frame(x1=5200,x2=2000,x3=1100,x4=1000,x5=1300,x6=45000,x7=34000,x8=115.0,x9=3.8)
> predict(dy,newdata,interval="confidence")
fit lwr upr
1 13718.67 13468.98 13968.36
>

> h=ts(read.csv("F://3.csv",header=TRUE))
> h
Time Series:
Start = 1
End = 56
Frequency = 1
X78
[1,] -58
[2,] 53
[3,] -63
[4,] 13
[5,] -6
[6,] -16
[7,] -14
[8,] 3
[9,] -74
[10,] 89
[11,] -48
[12,] -14
[13,] 32
[14,] 56
[15,] -86
[16,] -66
[17,] 50
[18,] 26
[19,] 59
[20,] -47
[21,] -83
[22,] 2
[23,] -1
[24,] 124
[25,] -106
[26,] 113
[27,] -76
[28,] -47
[29,] -32
[30,] 39
[31,] -30
[32,] 6
[33,] -73
[34,] 18
[35,] 2
[36,] -24
[37,] 23
[38,] -38
[39,] 91
[40,] -56
[41,] -58
[42,] 1
[43,] 14
[44,] -4
[45,] 77
[46,] -127
[47,] 97
[48,] 10
[49,] -28
[50,] -17
[51,] 23
[52,] -2
[53,] 48
[54,] -131
[55,] 65
[56,] -17
> plot(h,type="o")

> local({pkg <- select.list(sort(.packages(all.available = TRUE)),graphics=TRUE)
+ if(nchar(pkg)) library(pkg, character.only=TRUE)})
Warning message:
程辑包‘urca’是用R版本3.4.4 来建造的
> adf=ur.df(as.vector(h),type=c("drift"),selectlags=c("AIC"))
> summary(adf)

###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################

Test regression drift

Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)

Residuals:
Min 1Q Median 3Q Max
-96.191 -23.390 -0.581 18.446 133.241

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -9.4381 7.0489 -1.339 0.187
z.lag.1 -1.7837 0.2386 -7.476 9.65e-10 ***
z.diff.lag 0.1956 0.1379 1.418 0.162
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 50.89 on 51 degrees of freedom
Multiple R-squared: 0.7589, Adjusted R-squared: 0.7494
F-statistic: 80.25 on 2 and 51 DF, p-value: < 2.2e-16

Value of test-statistic is: -7.4761 27.9471

Critical values for test statistics:
1pct 5pct 10pct
tau2 -3.51 -2.89 -2.58
phi1 6.70 4.71 3.86

> acf(h)

> pacf(h)

> ar=sarima(h,1,0,4,details=F)
> ar
$fit

Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
Q), period = S), xreg = xmean, include.mean = FALSE, optim.control = list(trace = trc,
REPORT = 1, reltol = tol))

Coefficients:
ar1 ma1 ma2 ma3 ma4 xmean
-0.0957 -0.7605 -0.051 -0.2591 0.0706 -5.0886
s.e. 0.7318 0.7244 0.637 0.2013 0.1939 0.4252

sigma^2 estimated as 1850: log likelihood = -291.97, aic = 597.95

$degrees_of_freedom
[1] 50

$ttable
Estimate SE t.value p.value
ar1 -0.0957 0.7318 -0.1308 0.8965
ma1 -0.7605 0.7244 -1.0498 0.2988
ma2 -0.0510 0.6370 -0.0800 0.9365
ma3 -0.2591 0.2013 -1.2875 0.2038
ma4 0.0706 0.1939 0.3641 0.7173
xmean -5.0886 0.4252 -11.9668 0.0000

$AIC
[1] 8.73734

$AICc
[1] 8.814721

$BIC
[1] 7.954342

> ma=sarima(h,0,1,1,details=F)
> ma
$fit

Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1,
reltol = tol))

Coefficients:
ma1 constant
-1.0000 0.1275
s.e. 0.0452 0.4833

sigma^2 estimated as 3412: log likelihood = -303.77, aic = 613.53

$degrees_of_freedom
[1] 53

$ttable
Estimate SE t.value p.value
ma1 -1.0000 0.0452 -22.1390 0.000
constant 0.1275 0.4833 0.2638 0.793

$AIC
[1] 9.206399

$AICc
[1] 9.250355

$BIC
[1] 8.278733
> arma=sarima(h,1,1,1,details=F)
> arma
$fit

Call:
stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
Q), period = S), xreg = constant, optim.control = list(trace = trc, REPORT = 1,
reltol = tol))

Coefficients:
ar1 ma1 constant
-0.4893 -1.0000 0.1052
s.e. 0.1161 0.0469 0.2858

sigma^2 estimated as 2548: log likelihood = -296.27, aic = 600.53

$degrees_of_freedom
[1] 52

$ttable
Estimate SE t.value p.value
ar1 -0.4893 0.1161 -4.2127 0.0001
ma1 -1.0000 0.0469 -21.3207 0.0000
constant 0.1052 0.2858 0.3680 0.7143

$AIC
[1] 8.950118

$AICc
[1] 8.999838

$BIC
[1] 8.058619

> res=residuals(ar$fit)
> Box.test(res)

Box-Pierce test

data: res
X-squared = 0.0040697, df = 1, p-value = 0.9491

> plot(res*res)

> res<-residuals(ma$fit)
> res
Time Series:
Start = 1
End = 56
Frequency = 1
[1] -5.812742e-02 7.839872e+01 -4.955419e+01 3.066745e+01 6.646768e+00 -3.818017e+00 -1.493191e+00 1.448967e+01 -5.993782e+01 1.009045e+02 -3.947443e+01 -3.604570e+00 4.075719e+01 6.073798e+01 -8.076426e+01 -5.630655e+01
[17] 5.952211e+01 3.267028e+01 6.289877e+01 -4.376929e+01 -7.688972e+01 9.609734e+00 6.123687e+00 1.281064e+02 -1.026027e+02 1.160447e+02 -7.392804e+01 -4.288658e+01 -2.676745e+01 4.382151e+01 -2.561905e+01 1.050204e+01
[33] -6.774055e+01 2.380823e+01 7.222574e+00 -1.874297e+01 2.800543e+01 -3.305934e+01 9.500912e+01 -5.267336e+01 -5.347395e+01 5.982226e+00 1.856342e+01 2.163113e-01 8.018037e+01 -1.234786e+02 1.006553e+02 1.232087e+01
[49] -2.566963e+01 -1.438774e+01 2.537689e+01 -6.110995e-04 4.939921e+01 -1.289854e+02 6.746525e+01 -1.514136e+01

> Box.test(res)#

Box-Pierce test

data: res
X-squared = 13.335, df = 1, p-value = 0.0002606

> yc=sarima.for(h,10,1,1,1)

> yc$pred
Time Series:
Start = 57
End = 66
Frequency = 1
[1] 5.106162 -5.553160 -0.181147 -2.652867 -1.286844 -1.798532 -1.391503 -1.433980 -1.256525 -1.186677