Convert R code into matlab

  • #1
12
0

Main Question or Discussion Point

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:

Answers and Replies

  • #2
berkeman
Mentor
56,895
6,864
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.
 
  • #3
12
0
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.
 
  • #4
berkeman
Mentor
56,895
6,864
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?
 
  • #5
12
0
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.
 
  • #6
33,303
4,998
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)
 
  • #7
FactChecker
Science Advisor
Gold Member
5,391
1,957
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.
 
  • #8
DrDu
Science Advisor
6,023
755
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.
 
  • #9
DrDu
Science Advisor
6,023
755
Code:
require(tseries)
require(fracdiff)
require(matrixStats)
Yes, they can be found on the CRAN archive, including their descriptions.
 
  • #10

Related Threads on Convert R code into matlab

Replies
2
Views
2K
Replies
3
Views
1K
Replies
17
Views
2K
  • Last Post
Replies
1
Views
2K
  • Last Post
Replies
2
Views
1K
  • Last Post
Replies
2
Views
3K
  • Last Post
Replies
0
Views
1K
Replies
1
Views
6K
Replies
1
Views
2K
  • Last Post
Replies
1
Views
1K
Top