R Programming # 8 and 9

# 1. Assume the distribution of at as N (0, σa ) in the M A(1) model :

# Rt = at + 0.2at−1 , #σa = 0.025 #Generate 100 random returns (Rt )

# based on the above model. Plot the sample auto-

# correlation in different lags. Plot the actual autocorrelation with

# respect to different

# lags.

ma=function()

{

library(graphics);

library(MASS);

library(fBasics);

layout(matrix(c(1,2,3,4),2,2)) # optional 4 graphs/page

at=rnorm(1002, mean = 0, sd = 0.025);

i=0;

Rt=0;

for(i in 2:1000)

{

Rt[i]=at[i]+0.2*at[i-1];

}

#sample acf plot

lacfPlot(Rt, n = 5, lag.max = 25, type = c("returns", "values"), labels = TRUE)

lacfPlot(Rt, n = 10, lag.max = 25, type = c("returns", "values"), labels = TRUE)

lacfPlot(Rt, n = 15, lag.max = 25, type = c("returns", "values"), labels = TRUE)

lacfPlot(Rt, n = 20, lag.max = 25, type = c("returns", "values"), labels = TRUE)

lacfPlot(at, n = 5, lag.max = 25, type = c("returns", "values"), labels = TRUE)

lacfPlot(at, n = 10, lag.max = 25, type = c("returns", "values"), labels = TRUE)

lacfPlot(at, n = 15, lag.max = 25, type = c("returns", "values"), labels = TRUE)

lacfPlot(at, n = 20, lag.max = 25, type = c("returns", "values"), labels = TRUE)

}

 ===========================================

# 2. Plot the pdf of GEV when its parameter ξ = 0, 1 and -1. Superimpose

# the three graphs.

gev=function()

{

library(graphics);

library(MASS);

library(fBasics);

library(fExtremes);

layout(matrix(c(1,2,3,4),2,2)) # optional 4 graphs/page

x=seq(-10,10,length=10000);

d=dgev(x, xi = 0, mu = 0, beta = 1, log = FALSE)

plot(d, type = "l", col = "red", main = "Weibull Series with parameter = 0")

grid()

d=dgev(x, xi = 1, mu = 0, beta = 1, log = FALSE)

plot(d, type = "l", col = "blue", main = "Weibull Series with parameter = 1")

grid()

d=dgev(x, xi = -1, mu = 0, beta = 1, log = FALSE)

plot(d, type = "l", col = "green", main = "Weibull Series with parameter = -1")

grid()

d=dgev(x, xi = 0, mu = 0, beta = 1, log = FALSE)

plot(d, type = "l", col = "red", main = "")

grid()

par(new=TRUE)

d=dgev(x, xi = 1, mu = 0, beta = 1, log = FALSE)

plot(d, type = "l", col = "blue", main = "")

grid()

par(new=TRUE)

d=dgev(x, xi = -1, mu = 0, beta = 1, log = FALSE)

plot(d, type = "l", col = "green", main = "superimposition of parameter = 0(RED), 1(BLUE), -1(GREEN)")

grid()

}

# 3. Draw the mean excess plot for normal, exponential and pareto

# distribution.

mep=function()

{

library(graphics);

library(MASS);

library(fBasics);

library(fExtremes);

library(evir);

library(VGAM);

layout(matrix(c(1,2,3,4),2,2)) # optional 4 graphs/page

x=seq(-5,5,length=1000);

n=dnorm(x, mean = 0, sd = 1, log = FALSE);

e=dexp(x, rate=1);

p=dpareto(x, 2, 2, log=FALSE)

ln=dlnorm(x, meanlog = 0, sdlog = 1, log = FALSE)

meplot(n, omit = 3, labels = TRUE, main="MEP of normal distribution");

meplot(e, omit = 3, labels = TRUE, main="MEP of exponential distribution");

meplot(p, omit = 3, labels = TRUE, main="MEP of pareto distribution");

meplot(ln, omit = 3, labels = TRUE, main="MEP of lognormal distribution");

}

# 4. Consider the daily stock return of the Citigroup (tick symbol C)

# and the Standard and Poor’s 500 Composite index from January 2001 to

# December 2008. The data are simple returns and in the file

# d-csp0108.txt (three columns with date, C-rtn, SP-rtn).

# (a) Fit the generalized extreme value distribution to negative S&P

# log-returns, in percentages, with subperiods of 30 trading days.

# Write down the parameter estimates and their standard errors.

# Obtain a scatterplot and a QQ-plot of the residuals.

gevd=function()

{

library(graphics);

library(MASS);

library(fBasics);

library(fExtremes);

library(evir);

library(VGAM);

library(ismev);

library(dse1);

layout(matrix(c(1,2,3,4),2,2)) # optional 4 graphs/page

data=read.table("d-csp0108.txt",header=TRUE);

c=data[2];

sp=data[3];

c=cbind(c);

sp=cbind(sp);

#gev.fit(c)

#gev.fit(sp)

x=seq(-1,1,length=100);

#ft=gev(x, block=30)

r=residualStats(matrix(rnorm(200), 100,2), NULL)

print(r)

print(length(r))

x=r[1];

y=r[2]

plot(x,y)

qplot(r)

}

===================================================

pca=function()

{

data=read.table("m-pca5c-9003.txt",header=TRUE);

ibm=data[1];

hp=data[2];

intel=data[3];

ml=data[4];

msdw=data[5];

# Pricipal Components Analysis

# entering raw data and extracting PCs

# from the correlation matrix

fit <- princomp(data, cor=TRUE)

summary(fit) # print variance accounted for

loadings(fit) # pc loadings

plot(fit,type="lines") # scree plot

fit$scores # the principal components

biplot(fit)

# Varimax Rotated Principal Components

# retaining 5 components

library(psych)

fit <- principal(data, nfactors=5, rotate="cluster")

fit # print results

summary(fit) # print variance accounted for

loadings(fit) # pc loadings

plot(fit,type="lines") # scree plot

#fit$scores # the principal components

#biplot(fit)

}