# Convert R code into matlab

• MATLAB

## 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:

Related MATLAB, Maple, Mathematica, LaTeX News on Phys.org
berkeman
Mentor
I added code tags for readability -- please always use code tags to post programs here. Thanks.

So this is for schoolwork? Can you show us your Attempt at the Solution? We don't do your schoolwork for you here.

I added code tags for readability -- please always use code tags to post programs here. Thanks.

So this is for schoolwork? Can you show us your Attempt at the Solution? We don't do your schoolwork for you here.
It's not a schoolwork.
I would like to perform an anlysis of two time series x and y, my aim is the perform a particular regression approach, on them.

berkeman
Mentor
Okay, let me move it to the Programming forum then...

ADD -- Have you started the MATLAB code? Can you post what you have done so far?

Okay, let me move it to the Programming forum then...

ADD -- Have you started the MATLAB code? Can you post what you have done so far?
Not yet, i didn't understant what this codes do, i don't understand R programming language.

Mark44
Mentor
Not yet, i didn't understant what this codes do, i don't understand R programming language.
If you don't understand the R code, you will not be able to convert it to matlab code. If I were you I would start by looking at the R documentation.

The three lines below seem to be libraries that the code uses. If so, it would be useful to look at these libraries to see what the functions in them do, what parameters are required, and what kinds of values they return.
Code:
require(tseries)
require(fracdiff)
require(matrixStats)
It would also be useful to understand what this notation means. I don't know R, so I can't offer any advice, but I would be looking to see what DCCA_beta_avg is, and what the notation <- means.
Code:
DCCA_beta_avg<-function(y,x,smin,smax,step)
From the line below, it's possible that R uses <- for assignment operations, but I can't say for sure.
Code:
XTX<-var(x)*(length(x)-1)

FactChecker
Gold Member
Step 1 is to find the person who wrote the R code without a single comment and poke him in the eye. Do you even know what the input and the output looks like? If so, you may be able to write MATLAB code without knowing every detail of the R code.

DrDu
Why do you want to convert it to matlab? It is possible to run matlab from inside R and use the results for computation. This would be probably easier.

DrDu
Code:
require(tseries)
require(fracdiff)
require(matrixStats)
Yes, they can be found on the CRAN archive, including their descriptions.

robphy