Good morning everyone,
I don't want to bother you. I'm new at using R. :)
1. I was wondering if someone could help me figure out why I can't generate
the code to get a hyperbolic function.
2. My second question is, I generated the code. I don't have any problem
with other distributions but I still can't get the graphics displayed.
Here are the instructions for my exercise and here is the code I used:
**Instructions**
Project: hereafter the series of financial returns will be refered to as yt
and the series of fundamentals as xt. Here are the questions you need to
raise and answer:
Part 1: maximum likelihood estimation, student test and goodness of fit.
1. Consider the following model:
yt =μ+σεt, with εt a disturbance such that E[εt] = 0 and V[εt] = 1.
*Estimate the parameters of the following distributions by maximum
likelihood using the Yt:*
• Gaussian distribution N(0,1)
• Student-t distribution with parameter ν
• A mixture of Gaussian distribution MN(φ, μ1, σ1, μ2, σ2)
• A generalized hyperbolic distribution GH(λ, α, β, δ, μ).
2. *Test the parameters for their significance using a Student test. Are
all the parameters statistically significant?*
**My code:**
For the first three distributions, I answered well for questions one and
two which are in italics. But I have a problem with the last one.
library(readxl)
Data <- read_excel("Data_project.xlsx", sheet = 2, skip = 5, col_names =
FALSE)
#y variable to explain, Nikkei 225 index
y <- Data[16]
# x variable explanatory variable, Australian Dollar vs. US Dollar
x <- Data[30]
returns_y = diff(as.matrix(log(y)),1)
returns_x = diff(as.matrix(log(x)),1)
returnsy_std = scale(returns_y)
#1. Estimate the parameters by maximum likelihood
#Gaussian distribution N(0, 1)
mu=0
sigma=0.1
para0 = c(mu,sigma)
loglikG<-function(para,returns_y)
{
mu = para[1]
sigma = para[2]
print(para)
logl=sum(log(dnorm((returns_y-mu)/sigma, 0, 1)/sigma))
return(-logl)
}
loglikG(para0,returns_y)
fit <- optim(para0, loglikG, gr=NULL, returns_y, method="BFGS",hessian=T)
paraopt<- fit[["par"]]
Hessian = fit$hessian
invh = solve(Hessian)
t1= sqrt(invh[1,1])
t2=sqrt(invh[2,2])
testzmu = paraopt[1]/t1
testzvar = paraopt[2]/t2
print(testzmu)
print(testzvar)
#T-student
para_t=c(0,0.012,5)
loglikt <- function(para,returns_y){
mut=para[1]
sigmat=para[2]
nu=para[3]
m=-sum(log(dt((returns_y-mut)/sigmat, df=nu)/sigmat))
return(m)
}
loglikt(para_t,returns_y)
output_t= optim(para_t, loglikt, gr=NULL, returns_y, method="BFGS",
hessian=TRUE)
paraopt_t <- output_t[["par"]]
Hessian_t = output_t$hessian
invh_t = solve(Hessian_t)
t1_t= sqrt(invh_t[1,1])
t2_t=sqrt(invh_t[2,2])
t3_t= sqrt(invh_t[3,3])
testzmu_t = paraopt_t[1]/t1_t
testzvar_t = paraopt_t[2]/t2_t
testznu_t = paraopt_t[3]/t3_t
print(testzmu_t)
print(testzvar_t)
print(testznu_t)
#Mixture of Gaussian finding initial values
library(LaplacesDemon)
eps = 0.001
tolerance = 0.95
paraMG = c(-0.02,0.03,0.6,0.8,0.7)
likehoodMG <- function(para,returnsy_std)
{
muM12 = para[1:2]
sigmaM12 = para[3:4]
phi = para[5]
p = c(phi,1-phi)
LM=sum(log(dnormm(returnsy_std,muM12,sigmaM12,p=p)))
mean_w = p[1]*muM12[1] + p[2]*muM12[2]
var_w = p[1]*sqrt(sigmaM12[1]) + p[2]*sqrt(sigmaM12[2])
if((abs(mean_w)>tolerance) || (abs(var_w-1)>tolerance)){
return(NaN)
}
return(-LM)
}
likehoodMG(paraMG,returnsy_std)
outputMG = optim(paraMG, likehoodMG, gr=NULL, returnsy_std, method =
"L-BFGS-B", hessian=TRUE,
lower = c(eps,eps,-Inf,-Inf,eps), upper =
c(Inf,Inf,Inf,Inf,1-eps))
paraoptMG = outputMG[["par"]]
#Mixture of Gaussian
#(0.000345,0.023306)
paraM
=c(0.000345,0.023306,0.0253514,0.0010000,0.8715856,2.7857329,0.9659020)
likehoodM <- function(para,returns_y)
{
muM1 = para[1]
sigmaM = para[2]
muM12 = para[3:4]
sigmaM12 = para[5:6]
phi = para[7]
p = c(phi,1-phi)
LM=sum(log(dnormm((returns_y-muM1)/sigmaM,muM12,sigmaM12,p=p)/sigmaM))
return(-LM)
}
likehoodM(paraM,returns_y)
outputM = optim(paraM, likehoodM, gr=NULL, returns_y, method = "L-BFGS-B",
hessian=TRUE,
lower = c(-Inf,eps,eps,eps,-Inf,-Inf,eps), upper =
c(Inf,Inf,Inf,Inf,Inf,Inf,1-eps))
paraoptM = outputM[["par"]]
HM = outputM[["hessian"]]
invHM = solve(HM)
tm1 = sqrt(invHM[1,1])
tm2 = sqrt(invHM[2,2])
tm3 = sqrt(invHM[3,3])
tm4 = sqrt(invHM[4,4])
tm5 = sqrt(invHM[5,5])
tm6 = sqrt(invHM[6,6])
tm7 = sqrt(invHM[7,7])
testtmum = (paraoptM[1]-0)/tm1
testtsigmam = paraoptM[2]/tm2
testtmum1 = (paraoptM[3])/tm3
testtmum2 = paraoptM[4]/tm4
testtvarm1 = paraoptM[5]/tm5
testtvarm2 = paraoptM[6]/tm6
testtphi = paraoptM[7]/tm7
print(testtmum)
print(testtsigmam)
print(testtmum1)
print(testtmum2)
print(testtvarm1)
print(testtvarm2)
print(testtphi)
#A generalized hyperbolic distribution GH(??, ??, ??, ??, ?).
library(readxl)
Data <- read_excel("Data_project.xlsx", sheet = 2, skip = 5, col_names =
FALSE)
#y variable to explain, Nikkei 225 index
y <- Data[16]
# x variable explanatory variable, Australian Dollar vs. US Dollar
x <- Data[30]
returns_y = diff(as.matrix(log(y)),1)
returns_x = diff(as.matrix(log(x)),1)
returnsy_std = scale(returns_y)
#A generalized hyperbolic distribution GH(??, ??, ??, ??, ?).
library(ghyp)
para_gh= c(-0.00002,0.005,1,0,1,0.1,0.1) # mu,delta,alpha,beta,lambda
loglikGH <- function(para,returns_y){
mugh=para[1]
sigmagh=para[2]
alpha=para[3]
beta = para[4]
delta=para[5]
chi=para[6]
lamda=para[7]
if(delta < abs(chi)){
return(10000)
}else{
return(-sum(log(dghyp(((returns_y-mugh)/sigmagh),object =
ghyp(alpha,beta,delta,chi,lamda))/sigmagh)))
}
}
loglikGH(para_gh,returns_y)
eps = 0.001
outputGH= optim(para_gh, loglikGH, gr=NULL, returns_y, method = "L-BFGS-B",
hessian=TRUE,
lower = c(-Inf,eps,eps,eps,-Inf,-Inf,eps), upper =
c(Inf,Inf,Inf,Inf,Inf,Inf,Inf))
paraoptGH = outputGH[["par"]]
HGH = outputGH[["hessian"]]
invHGM = solve(HGH)
tm1_H = sqrt(invHGM[1,1])
tm2_H = sqrt(invHGM[2,2])
tm3_H = sqrt(invHGM[3,3])
tm4_H = sqrt(invHGM[4,4])
tm5_H = sqrt(invHGM[5,5])
tm6_H = sqrt(invHGM[6,6])
tm7_H = sqrt(invHGM[7,7])
testtmuH = (paraoptGH[1]-0)/tm1_H
testtsigmaH = paraoptGH[2]/tm2_H
testtalpha = (paraoptGH[3])/tm3_H
testtbetha = paraoptGH[4]/tm4_H
testtdelta = paraoptGH[5]/tm5_H
testtvchi = paraoptGH[6]/tm6_H
testtlambda = paraoptGH[7]/tm7
print(testtmuH)
print(testtsigmaH)
testtalpha
testtbetha
testtdelta
testtvchi
testtlambda
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
