典型例题
HW2——5
HW3——3，4
HW4——
EXAM14,15——VAR
藏族美食
library(quantmod)
getSymbols("AAPL",from="2008-01-03",to="2015-01-28")
dim(AAPL)
head(AAPL)
tail(AAPL)
chartSeries(AAPL,theme="white") % Obtain time plot of closing price and trading volume
da <- read.csv("d-vix0411.csv",header=T)%有时候不用header参数，看head(da)的数据形式
da <- read.table("",header=T)
=diff(log(AAPL$AAPL.Adjusted)) % Compute log returns
unrate <- as.numeric(UNRATE[,1]) % U a regular vector, instead of an“xts”object
library(fBasics)
ts.plot(ibm,main="Monthly IBM simple returns: 1968-2015") % Time plot
basicStats(ibm)
apply(rtn,2,basicStats) ## This command says apply "basicStats" to every 信封封面设计
columns in "rtn"
d4=density(ibm)
plot(d4$x,d4$y,type='l',xlab="rtn",ylab='AAPL') ##density图
t.test(lnIBM) %% Test mean=0 vs mean .not. Zero也可以用basicStats命令去 计算
normalTest(lnIBM,method=’jb’)高贵近义词 %testing normality of financial return ries
s3=skewness(lnIBM); T <- length(lnIBM)
tst <- s3/sqrt(6/T)赢在执行 % test skewness
pv <- 2*pnorm(tst) %calculate p value
k4 <- kurtosis(lnIBM)
tst <- k4/sqrt(24/T) % test excess kurtosis
mu <- mean(sbux); v1 <- var(sbux) %%prediction 或者可以用t-test的信息
lcl <- mu-1.96*sqrt(v1)
ucl <- mu+1.96*sqrt(v1)
c(lcl,ucl)
correlation
cor(sp,ibm)
[1] 0.5785249
cor(sp,ibm,method="kendall")月加青念什么 %random copy
[1] 0.4172056
cor(sp,ibm,method="spearman") %rank correlation
[1] 0.58267
cor(rank(ibm),rank(sp))
[1] 0.58267
x[1,] 孕妇梦见甘蔗 %show the first row of the data
X[,4:5] %show the 4and5 columns of the data
y=ts(x[,3],frequency=252,start=c(2004,1)) <== Create a time-ries object in R
plot(y,type='l',xlab='year',ylab='rtn')
par(mfcol=c(2,1)) <== To put two plots on a single page
hist(y,nclass=50)%直方图
tdx=(c(1:615)+11)/12+1959 %创建time plot 的横坐标
plot(tdx,xt,xlab='year',ylab='temp',type='l')
plot(tdx[-1],zt,xlab='year',ylab='diff(temp)',type='l') ##注意difference使得横坐标年份变化
int=cbind(x[300:914,4],y[,4]) <== Line up the two TB rates
acf(ibm)
芭比兔m1 <- acf(ibm)
names(m1)
m1$acf
m2 <- pacf(ibm)
st(ibm,lag=10) % Box-Pierce Q(m) test，序列是否自相关
st(ibm,lag=10,type='Ljung') % Ljung-Box Q(m) test
m1=lm(r3~r1) % Fit a regression model with likelihood method
● AR,MA,ARMA,ARIMA model
m1=ar(x,method='mle') % Automatic AR fitting using AIC criterion
names(m1)
m1$order
plot(m1$resid,type='l')
或者acf(m1$residuals) % Plot residuals of the fitted model
st(m1$resid,lag=10,type='Ljung') % Model checking
m2=arima(x,order=c(3,0,0)) % Another approach with order given
tsdiag(m2) % obtain 3 plots of model checking
tsdiag(m2,gof=24) %Model checking with 24 lags of residuals ACF
c1 <- c(NA,NA,NA,0,0,NA,0,0,0,NA,NA,NA) %Fixing coefficients to zero.根据t值将不显著的参数设置为0，注意一旦看出t值不显著，就要做这一步
m3 <- arima(unrate,order=c(11,0,0),fixed=c1)
jan <- rep(c(1,rep(0,11)),54) %%定义虚拟变量——联系下面handle outliers
tsdiag(m3,gof=24)
source（backtest.R）
backtest(m2,lhp,200,1)审核凭证%比较两个model哪个更好
backtest(n5,dB,orig=550,xre=an=F)
Dealing with the constant term. If there is any differencing, no constant is ud.
The subcommand an=F in the arima command.
有an=F的时候估计的模型里面没有常数项
% Further analysis of the fitted model.
p1=c(1,-m2$coef[1:3])
roots=polyroot(p1)
roots(大于一)
Mod(roots)
predict(m2,8) % Prediction 1-step to 8-step ahead and intervals
ucl <- p3$pred+1.96*p3$ %%interval
lcl <- p3$pred-1.96*p3$
cint <- cbind(lcl,ucl)
print(cint)
### Handling outliers
> which.min(m3$residuals) ### Locating the largest outlier
[1] 23
> I23 <- rep(0,819) ## Create indicator variable
> I23[23]=1
> c2 <- c(c1,NA) ## Include the new parameter
> m4 <- arima(unrate,order=c(11,0,0),fixed=c2,xreg=I23) %%有时候xreg还有其他X变量
> tsdiag(m4,gof=24)
library(fUnitRoots) %unit-root test
adfTest(gdp,lag=4,type=c("c")) #用ar(diff(gdp),method=”mle”)命令去获取lag
library(forecast)
Auto.arima(x) ##用此命令之后要看acf图改进模型，再根据t值改进
如果说auto出来的有with zero mean，一定要加上an=F
##Seasonal Time Series
> x=ts(scan(""),frequency=4,start=c(1960,1)) % create a time ries object
> plot(x) % Plot data with calendar time
> y=log(x) % Natural log transformation
> plot(y) % plot data 看是否有季节因素影响趋势
> c1=paste(c(1:4))
