- #1

Tilfani

- 12

- 0

## Homework Statement

I would like to convert this R code into matlab. This code perform a calculation of beta coefficient based on Detrended cross correlation analysis.

## Homework Equations

Code:

```
require(tseries)
require(fracdiff)
require(matrixStats)DCCA_beta_avg<-function(y,x,smin,smax,step){
XTX<-var(x)*(length(x)-1)
betas<-rep(0,(smax-smin)/step+1)
for(s in seq(smin,smax,by=step)){
betas[(s-smin)/step+1]<-DCCA_beta_sides(y,x,s)
}
DCCA_beta<-mean(betas)
DCCA_res<-(y-DCCA_beta*x)-mean(y-DCCA_beta*x)
DCCA_sigma2<-sum(DCCA_res^2)/(length(DCCA_res)-2)
DCCA_SE<-sqrt(DCCA_sigma2/XTX)
DCCA_R2<-1-var(DCCA_res)/var(y)
OLS_beta<-lm(y~x)$coefficients[2]
OLS_res<-(y-OLS_beta*x)-mean(y-OLS_beta*x)
OLS_sigma2<-sum(OLS_res^2)/(length(OLS_res)-2)
OLS_SE<-sqrt(OLS_sigma2/XTX)
OLS_R2<-1-var(OLS_res)/var(y)
return(c(OLS_beta,OLS_SE,OLS_R2,DCCA_beta,DCCA_SE,DCCA_R2))
}DCCA_beta<-DCCdccafunction(y,x,s){
xx<-cumsum(x-mean(x))
yy<-cumsum(y-mean(y))
t<-1:length(xx)
F2sj_xy<-runif(floor(length(xx)/s))
F2sj_xx<-F2sj_xy
for(ss in seq(1,(floor(length(xx)/s)*s),by=s)){
F2sj_xy[(ss-1)/s+1]<-sum((summary(lm(xx[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(yy[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)
F2sj_xx[(ss-1)/s+1]<-sum((summary(lm(xx[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(xx[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)
}
beta<-mean(F2sj_xy)/mean(F2sj_xx)
return(beta)
}DCCA_beta_F<-function(y,x,s){
xx<-cumsum(x-mean(x))
yy<-cumsum(y-mean(y))
t<-1:length(xx)
F2sj_xy<-runif(floor(length(xx)/s))
F2sj_xx<-F2sj_xy
F2sj_yy<-F2sj_xy
for(ss in seq(1,(floor(length(xx)/s)*s),by=s)){
F2sj_xy[(ss-1)/s+1]<-sum((summary(lm(xx[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(yy[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)
F2sj_xx[(ss-1)/s+1]<-sum((summary(lm(xx[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(xx[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)
F2sj_yy[(ss-1)/s+1]<-sum((summary(lm(yy[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(yy[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)
}
beta<-mean(F2sj_xy)/mean(F2sj_xx)
return(c(beta,mean(F2sj_xx),mean(F2sj_yy)))
#return(c(beta,sum(F2sj_xx),sum(F2sj_yy)))
}DCCA_beta_SE<-function(y,x,s){
r<-DCCA_beta_F(y,x,s)
beta<-r[1]
yhat<-beta*x
alpha<-mean(y)-beta*mean(x)
res<-y-yhat
residuals<-res-mean(res)
resres<-cumsum(residuals-mean(residuals))
F2sj_res<-runif(floor(length(residuals)/s))
t<-1:length(resres)
for(ss in seq(1,(floor(length(residuals)/s)*s),by=s)){
F2sj_res[(ss-1)/s+1]<-sum((summary(lm(resres[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(resres[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)
}
#SE<-mean(residuals^2)/((length(residuals)-2)*r[2])
SE<-mean(F2sj_res)/((length(residuals)-2)*r[2])
SE_a<-(mean(F2sj_res)/r[2])*(sum(x^2)/(length(residuals)*(length(residuals)-2)))
R<-1-mean(F2sj_res)/(r[3])
return(c(alpha,sqrt(SE_a),beta,sqrt(SE),R))
}DCCA_beta_SE_F<-function(y,x,s){
r<-DCCA_beta_F(y,x,s)
beta<-r[1]
yhat<-beta*x
alpha<-mean(y)-beta*mean(x)
res<-y-yhat
residuals<-res-mean(res)
res_R<-y-x
resres<-cumsum(residuals-mean(residuals))
resres_R<-cumsum(res_R)
F2sj_res<-runif(floor(length(residuals)/s))
F2sj_res_R<-runif(floor(length(res_R)/s))
t<-1:length(resres)
for(ss in seq(1,(floor(length(residuals)/s)*s),by=s)){
F2sj_res[(ss-1)/s+1]<-sum((summary(lm(resres[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(resres[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)
F2sj_res_R[(ss-1)/s+1]<-sum((summary(lm(resres_R[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(resres_R[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)
}
#SE<-mean(residuals^2)/((length(residuals)-2)*r[2])
#SE<-mean(F2sj_res)/((length(residuals)-2)*r[2])
#SE<-mean(F2sj_res)/((length(F2sj_res)-2)*r[2]) #controlling for uncertainty connected to scales (higher scales have higher uncertainty due to lower number of blocks)
SE<-mean(F2sj_res)/(ceiling(length(residuals)/s)*r[2]) #loosing d.f. due to fitting a and b in each box
#SE_a<-(mean(F2sj_res)/r[2])*(sum(x^2)/(length(residuals)*(length(residuals)-2)))
#SE_a<-(mean(F2sj_res)/r[2])*(sum(x^2)/(length(F2sj_res)*(length(F2sj_res)-2))) #controlling for uncertainty connected to scales (higher scales have higher uncertainty due to lower number of blocks)
SE_a<-(mean(F2sj_res)/r[2])*(sum(x^2)/(length(residuals)*ceiling(length(residuals)/s))) #loosing d.f. due to fitting a and b in each box
R<-1-mean(F2sj_res)/(r[3])
#SSR_U<-sum(residuals^2)
SSR_U<-sum(F2sj_res)
#SSR_R<-sum((y-x)^2) #specific null: alpha=0, beta=1
SSR_R<-sum(F2sj_res_R)
#F_stat<-((SSR_R-SSR_U)/(SSR_U))*((length(residuals)-2)/2)
#F_stat<-((SSR_R-SSR_U)/(SSR_U))*((length(F2sj_res)-2)/2) #controlling for uncertainty connected to scales (higher scales have higher uncertainty due to lower number of blocks)
F_stat<-((SSR_R-SSR_U)/(SSR_U))*(ceiling(length(residuals)/s)/2) #loosing d.f. due to fitting a and b in each box
F_p<-pf(F_stat,2,length(F2sj_res)-2,lower.tail=FALSE)
return(c(alpha,sqrt(SE_a),beta,sqrt(SE),R,F_stat,F_p))
}DCCA_beta_s<-function(y,x,smin,smax,step){
results<-matrix(rep(0,6*((smax-smin)/step+1)),ncol=6)
for(s in seq(smin,smax,by=step)){
beta<-DCCA_beta_SE(y,x,s)
results[((s-smin)/step+1),1]<-s
results[((s-smin)/step+1),2]<-beta[1]
results[((s-smin)/step+1),3]<-beta[2]
results[((s-smin)/step+1),4]<-beta[3]
results[((s-smin)/step+1),5]<-beta[4]
results[((s-smin)/step+1),6]<-beta[5]
}
return(results)
}DCCA_beta_s_F<-function(y,x,smin,smax,step){
results<-matrix(rep(0,10*((smax-smin)/step+2)),ncol=10)
for(s in seq(smin,smax,by=step)){
beta<-DCCA_beta_SE_F(y,x,s)
results[((s-smin)/step+1),1]<-s
results[((s-smin)/step+1),2]<-beta[1]
results[((s-smin)/step+1),3]<-beta[2]
results[((s-smin)/step+1),4]<-2*pnorm(abs(beta[1]/beta[2]),lower.tail=FALSE)#p-value for null=0
results[((s-smin)/step+1),5]<-beta[3]
results[((s-smin)/step+1),6]<-beta[4]
results[((s-smin)/step+1),7]<-2*pnorm(abs((beta[3]-1)/beta[4]),lower.tail=FALSE)#p-value for null=1
results[((s-smin)/step+1),8]<-beta[5]
results[((s-smin)/step+1),9]<-beta[6]
results[((s-smin)/step+1),10]<-beta[7]
}
#results[(smax-smin)/step+2,2]<-mean(results[1:(dim(results)[1]-1),2])#A
#results[(smax-smin)/step+2,5]<-mean(results[1:(dim(results)[1]-1),5])#B
results[(smax-smin)/step+2,2]<-sum(results[1:(dim(results)[1]-1),2]*results[1:(dim(results)[1]-1),8])/sum(results[1:(dim(results)[1]-1),8])#A as R2(s) weighted
results[(smax-smin)/step+2,5]<-sum(results[1:(dim(results)[1]-1),5]*results[1:(dim(results)[1]-1),8])/sum(results[1:(dim(results)[1]-1),8])#B as R2(s) weighted
results[(smax-smin)/step+2,3]<-sqrt((sum(x^2)/length(x))*sum((y-results[(smax-smin)/step+2,2]-results[(smax-smin)/step+2,5]*x)^2)/((length(y)-dim(results)[1]+1)*sum((x-mean(x))^2)))#SE_A
results[(smax-smin)/step+2,4]<-2*pnorm(abs(results[(smax-smin)/step+2,2]/results[(smax-smin)/step+2,3]),lower.tail=FALSE)#p-value for null=0
results[(smax-smin)/step+2,6]<-sqrt(sum((y-results[(smax-smin)/step+2,2]-results[(smax-smin)/step+2,5]*x)^2)/((length(y)-dim(results)[1]+1)*sum((x-mean(x))^2)))#SE_B
results[(smax-smin)/step+2,7]<-2*pnorm(abs((results[(smax-smin)/step+2,5]-1)/results[(smax-smin)/step+2,6]),lower.tail=FALSE)#p-value for null=1
results[(smax-smin)/step+2,8]<-1-sum((y-results[(smax-smin)/step+2,2]-results[(smax-smin)/step+2,5]*x)^2)/sum(y^2)#R2
results[(smax-smin)/step+2,9]<-((length(x)-2)/2)*((results[(smax-smin)/step+2,8]-(1-(sum((y-x)^2))/(sum((y-mean(y))^2))))/(1-results[(smax-smin)/step+2,8]))#F_test_R2_based
results[(smax-smin)/step+2,10]<-pf(results[(smax-smin)/step+2,9],2,length(y)-2,lower.tail=FALSE)#F_test p_val
return(results)
}
```

## The Attempt at a Solution

I would like to obtain a MATLAB code.[/B]

Last edited by a moderator: